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Chapter I 


Background and Overview of Document 


I . Introduction 

The following summarizes the progress toward completion of a 
comprehensive diagnostic objective analysis system based upon the 
calculus of variations. The approach was to first develop the 
objective analysis subject to the constraints that the final 
product satisfies the five basic primitive equations for a dry 
inviscid atmosphere: the two nonlinear horizontal momentum 
equations, the continuity equation, the hydrostatic equation, and 
the thermodynamic equation. Then, having derived the basic model, 
there would be added to it the equations for moist atmospheric 
processes and the radiative transfer equation. 

The remainder of this section is organized into subsections as 
follows. A brief review of the problem design and progress to the 
beginning of the current completed grant period is given in 1.1. 
This first period is designated as Phase I. Then subsection 1.2 
summarizes progress under this grant (Phase II) and references 
following chapters that give results in greater detail. 

The reader should pay particular attention to the findings 
presented in Chapter V. A conceptual error within one of the 
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mathematical derivations was discovered while preparing Chapter V. 
The outcome is that some of the negative conclusions regarding the 
meaningful incorporation of satellite thermodynamic data into the 
meteorological data mainstream are no longer valid. 

1.2. Background of Initial Project and Phase I Summary 

By 1981, most quantitative satellite data for the troposphere 
consisted of TIROS-N temperature retrievals (Smith and Woolf, 
1976) , some cloud wind vectors, early VAS soundings (Smith, 1970; 
Chesters, et al., 1981), and initial SASS surface winds (Jones, et 
al., 1979). Although these data revealed new and interesting 
phenomena, it was becoming apparent that satellite data could not 
be used to quantify the dynamics of the troposphere without being 
supplemented by conventional data. 

The tools available to dynamically assimilate satellite data 
with conventional data consisted mostly of initialization methods 
for numerical prediction models. For example, normal mode 
initialization (Baer, 1977; Machenhaur, 1977) improved the dynamic 
coupling between the observations and numerical models. Early 
comparisons between numerical predictions with and without 
satellite temperature data using initialization schemes of the time 
were inconclusive (Phillips, 1976; Tracton et al., 1981; Ghil et 
al., 1979). The impact of satellite data was highly dependent on 
the capability of the analysis and forecast system used for the 
impact testing. 
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There was, therefore, a need for other methods for blending 
satellite data with conventional data. The direct variational 
methods presented by Sasaki (1958, 1970) are applicable to a 

diagnostic data assimilation method. Achtemeier (1975) applied 
Sasaki's method to blend diverse conventional data to satisfy the 
primitive equations. This work revealed the potential of the 
direct variational method as well as some of the inherent 
difficulties. The purpose of the variational objective analysis 
project funded by NASA was to extend the Achtemeier method to blend 
satellite data with conventional data. 

The general goals of the project which began in 1982 (see 
Achtemeier, et al., 1986) were to: 

(1) Modify the Achtemeier (1975) numerical variational model for 
the assimilation or merging of weather data (constraints: the 
nonlinear set of dynamical equations for atmospheric flow plus 
the radiative transfer equation) collected remotely from space 
with conventional weather data. 

(2) Translate the variational theory into a practical model that 
can be used by other scientists within the NASA Global Scale 
Atmospheric Processes Research Program and elsewhere. 

(3) Design the variational assimilation to be independent from 
numerical forecast models - that is, make it diagnostic. This 
does not preclude the method from being used to develop 
initial fields for numerical models nor the blending of model 
forecast fields with satellite and conventional data. 

(4) Apply the variational model to estimate minimal data 
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requirements (data type, sampling frequency, density, 
accuracy) for accurate diagnosis of atmospheric structure that 
could be useful for the design of future satellite-borne 
instruments . 

The project began with a review of the Achtemeier (1972, 1975, 
1979) applications of the Sasaki (1958, 1970) variational methods 
by a mathematical consultant. This review confirmed that the 
applications of Sasaki's methods were mathematically sound. We 
then commenced the development of a general numerical variational 
model that includes the following dynamical constraints: the two 
nonlinear horizontal momentum equations, the hydrostatic equation, 
the continuity equation, the thermodynamic equation for moist 
convective processes, a moisture conservation equation, and the 
radiative transfer equation. From the experience gained from the 
previous studies, it was known that the application of the direct 
variational method to the above dynamic constraints presented three 
formidable mathematical problems. 

) 

(1) Courant (1936) showed that the number of subsidiary conditions 
(dynamic constraints) must be at least one less than the 
number of adjustable dependent variables. Inclusion of the 
same number of constraints as dependent variables, such as the 
five primitive equations with five dependent variables, 
overdetermines the problem and a solution is not guaranteed. 
Achtemeier (1975) originally attempted to circumvent this 
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problem through a parameterization of the tendency terms of 
temperature and the velocity components that required the 
exact solution of the integrated continuity equation. 
However, this method (a variational adjustment within a 
variational adjustment) was considered a failure after an 
extensive analysis (Achtemeier, 1979) found unrealistically 
large velocity component tendencies where actual velocity 
changes over a 12-hr period were small. 

(2) Application of the direct variational method to the local 
tendencies of temperature and the horizontal velocity 
components yield terms in the Euler-Lagrange equations that 
are local tendencies of Lagrange multipliers. Boundary 
conditions for these terms are unknown. The problem of time 
consistency in variational problems has been examined by Lewis 
(1980, 1982) and Lewis et al (1983) . More recently, the time 
consistency problem has been found more tractable through use 
of the adjoint method (Lewis and Derber, 1985; Talagrand and 
Courtier, 1987) . 

(3) The Euler-Lagrange equations produced by Sasaki's method 
include the dynamical constraints plus variational equations 
which are equal in number to the number of dependent variables 
to be blended and equal in complexity to the complexity of the 
original constraints. This set of complicated nonlinear 
partial differential equations presents formidable programming 
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problems and is difficult to solve by traditional methods. 
Achtemeier (1975) developed a cyclical solution method which, 
upon scaling the equations and expressing the terms in powers 
of the Rossby number, places higher order terms in forcing 
functions, solves the equations as a linear set and updates 
the higher order terms with the new solutions. The method, 
which converges rapidly, was modified for this project. 

To better assure progress toward the development of the 
general variational analysis model, the task was divided into four 
variational models of increasing complexity. MODEL I, including 
the two horizontal momentum equations, the hydrostatic equation, 
and an integrated continuity equation as dynamical constraints, 
addresses the problems that arise from applying the direct 
variational method to local tendency terms. The cyclical solution 
method developed by Achtemeier (1975) is applied to the MODEL I 
finite differencing scheme. A nonlinear vertical coordinate reduces 
the need to vertically interpolate satellite-derived temperatures. 
MODEL II, which includes MODEL I plus the thermodynamic equation 
for a dry atmosphere, focuses on the overdetermined solution 
problem posed by Courant. MODEL III includes MODEL II plus the 
radiative transfer equation as a dynamical constraint and the 
radiance as a dependent variable. MODEL IV incudes moisture and 
its parameter izations. 

The theoretical development and initial evaluations of MODEL 
I are given by Achtemeier (1986a), Achtemeier et al. (1986a, 
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1986b) , and Kidder and Achtemeier (1986) . In MODEL I, the tendency 
terms for the horizontal velocity components were separated into 
advective and developmental components; the advective component 
explains local changes caused by disturbances moving in a basic 
current and the developmental component explains local tendencies 
caused by changes in the structures or intensities of the 
disturbances. The advective components were incorporated into the 
dynamic equations through time to space conversion, and the 
developmental components became dependent variables. 

MODEL I successfully diagnosed local tendencies of the 
horizontal velocity components. Diagnosed local tendencies of the 
horizontal velocity components compared favorably with observed 
three hour tendencies, whereas local tendencies calculated from 
direct substitution of unadjusted fields of data into the dynamic 
constraints did not. 

1.2 Summary of Results under Phase II 

Theoretical development, programming, and evaluation of MODEL 
II was completed early in 1987. MODEL II included the five 
primitive equations for a dry atmosphere as dynamic constraints 
ie. , MODEL I plus the thermodynamic equation. The partition of the 
horizontal velocity tendencies with the definition of the 
developmental components as dependent variables increased the 
number of dependent variables from five to seven. Therefore, we 
have solved the problem of overdetermining the solution (Courant, 
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1936) . 

Evaluation of MODEL II revealed that the tendency term 
formulations remained stable (Achtemeier, 1988c) . However, first 
order terms that contain the divergence adjustment cancel out in 
the cyclical solution formulations. The divergence adjustment must 
then be carried in second order terms and through other variables. 
Formulating the continuity equation as an ancillary constraint to 
be satisfied at each cycle - a variational model within a 
variational model - and "nudge" the solution does not always 
provide the desired dynamic balance with the divergent part of the 
wind. Details of the MODEL II evaluation are presented in Chapter 
II. 

Further examination of MODEL II found that the divergent 
component of the wind can be made to satisfy the continuity 
equation subject to the satisfaction of an additional mathematical 
constraint. That constraint requires the adjusted wind field to 
satisfy a particular solution of the integrated vorticity equation. 
The particular solution satisfies the conditions that the sum of 
the terms of the vorticity equation vanish when integrated from the 
surface to the top of the domain and that the divergence vanishes 
when integrated from the surface to the top of the domain. This 
latter condition satisfies the integrated continuity equation. 

Both MODEL I and MODEL II were reformulated mathematically and 
reprogrammed into a new version MODEL IIB. MODEL IIB is the 
subject of Chapter III. It forms the new basic 5-primitive 
equation variational objective analysis model that must be 
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satisfied before any additional dynamic constraints can be added to 
increase the complexity of the objective analysis technique. The 
successful completion of MODEL IIB therefore, is the focus of the 
Phase II research effort. 

The need to rederive the basic variational models also changed 
the focus of this research. The original objective to develop an 
analysis tool for the incorporation of satellite data into the 
mainstream of meteorological data was changed from applied research 
to basic research. Rather than take an existing method and modify 
it to produce a product, it became necessary to prove that the 
variational method fundamentally was workable. 

Several other research efforts aimed at expanding the 
variational objective model beyond the basic model were pursued 
commensurate with the development of MODEL IIB. A variationally 
constrained temperature profile analysis was derived with the 
dynamic constraints being the four radiative transfer equations for 
the microwave channels used to sample within the troposphere. 
Rawinsonde temperature and the four brightness temperatures are the 
adjustable variables. Chapter IV presents the results from this 
study. Then, in Chapter V, the radiative transfer equations were 
combined with the equations for a geostrophic, hydrostatic 
atmosphere as part of a theoretical study to determine the impact 
of the radiative transfer equations upon the variational analysis. 

Additional effort was directed to the goal was to make the 
variational analysis more responsive to the original observations. 
The gridding of meteorological data is the first step in performing 
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objective analysis and is done independently from the variational 
adjustments. Therefore our first efforts in achieving this goal 
was to improve upon the objective analysis method. Chapter VI 
presents a study of the effects of influence radius upon an 
objective analysis and Chapter VII presents a modification of a 
successive corrections technique for improved derivative 


calculations. 
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1. Introduction 

The MODEL II variational data assimilation model is the second 
of the four variational models designed to blend diverse 
meteorological data into a dynamically constrained data set. MODEL 
II differs from the MODEL I developed during Phase I in that it 
includes the thermodynamic equation as the fifth dynamical 
constraint. 

Thus MODEL II includes all five of the primitive equations 
that govern atmospheric flow for a dry atmosphere. The reason for 
delaying the introduction of the thermodynamic equation until MODEL 
II is as follows. Courant (1936) showed that the number of 
subsidiary conditions (dynamic constraints) must be at least one 
less than the number of adjustable dependent variables. The five 
primitive equations form a closed set of equations with five 
dependent variables. Inclusion of the same number of constraints 
as dependent variables overdetermines the problem and a solution is 
not guaranteed. Achtemeier (1975) attempted to circumvent this 
problem through a parameterization of the tendency terms of the 
velocity components and the temperature that required the exact 
solution of the integrated continuity equation. This method, a 
variational adjustment within a variational adjustment, was 
considered a failure after an extensive analysis (Achtemeier, 1979) 
found unrealistically large velocity component tendencies where 
actual velocity changes over a 12-hr period were small. 

The approach taken in the development of MODEL I was to make 
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possible the inclusion of the five primitive equations by 
increasing the number of dependent variables. We defined two new 
dependent variables, the developmental components of the horizontal 
velocity tendencies, which increased the number of dependent 
variables from five to seven. Though this solves the problem of 
the number of subsidiary conditions, the extent of internal 
coupling among the variables and within the equations could not be 
determined fully until the development and evaluation of MODEL II. 

2. MODEL II: Thermodynamic Equation as a Dynamic Constraint 


Upon defluxing and omitting the dissipation term of the 
thermodynamic equation in Anthes and Warner (1978) , the 
thermodynamic equation as it appears as a dynamical constraint in 
MODEL II is, 


dT Mmf ,,dT^„dT,. dT RTa 

-5- + m{u-z- + v-r-) + 0-5 

at ox ay da 


( 1 ) 


c pP c p 

The omega-term (term 4) of the thermodynamic equation can be 
transformed into the nonlinear sigma coordinate system through the 
definition, 


» P Pu 
P*-Pu 


o-p (p-p*) 3 +o 


( 2 ) 
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where the superscript, *, and the subscript, u, identify, 
respectively, the variables at the reference pressure level and at 
the top of the model atmosphere. For more information on the 
nonlinear vertical coordinate system, refer to Appendix A.l. 
Furthermore , 

p-d-o'-Effil] (p,-pT 3 (3) 

P -Pu 


where the subscript, s, refers to quantities measured at the 
surface. We differentiate (2) with respect to time. If 


and 


J~ [3p (p-p*) 2 +a] ( p-p *) 

then we may define two coefficients such that 


and 


<? 3 ' 


P-P 

JP 


* 


(5) 


( 6 ) 



f P-Pl 
K Ps~P* 


\4 

/ 


( 7 ) 
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for p>p*, and 


1 _ P*-Pu 
up a*p 


( 8 ) 


and 


<? 4 -0 ( 9 ) 

for p < p*. 

The thermodynamic equation in the nonlinear sigma coordinates 
is, upon substitution for the omega-term, 


dT w 

dt 


/ & T w 
+m(u-^+v 

ox 


87 V 

8y 


+d 


do 


-^(qr 3 d + g 4 to s )-J2--0 


( 10 ) 


Here the subscript, W, refers to the whole temperature, T y = T R + 
T, where T R is a reference temperature for the layer and is always 
in hydrostatic balance and T is the departure from the reference 
temperature that is subject to adjustment within the variational 
model. Substitution for the whole temperature yields the 
thermodynamic equation in the adjustable part of the temperature, 


dt ox oy do 


J£--£L( Tr -T) (q 3 a +q4 u> s ) + 6^-JL - o (11) 


dT 
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Now nondimensionalize the thermodynamic equation. Letting 

a 

u-^Uu 1 , v=Uv', At-d/OAt 7 , Ax-LAx 7 , 

T r -Qt' r - ( gH/R) t' r , A ( gH/R) ( F/R 0 ) bT 1 ( 12 ) 

p=Pp', 6~ (C/L) o' , (0 s ~ ( PC/ L) (o'g 

and dividing through by (C/L) (gH/R) (F/R 0 2 ) , the nondimensionalized 
thermodynamic equation with primes suppressed is, 


-R 


dT , dT dT, . dT, 
3- +m ( Uy- + v— ) +o 3- ] 
ofc ox dy do 


i? 


p 


(^t'r+T) (g 3 d+g 4 o s ) - [ 


LRRl 

CgHF 



( 13 ) 


Dividing by the additional R 0 renders (13) into the same order of 
magnitude as the other dynamic equations of MODEL II. In addition, 
it can be shown that the two terms that include T R combine to form 
the static stability, 



i 

do 



( 14 ) 


Therefore, the thermodynamic equation reduces to 


dT 


dT . „dT . , dT R 


R a [ + m ( u + 6 ( - -fL rg, ) ] 


dt 


dx dy 


do 


, LRR a n r R a 


( 15 ) 
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Next, the thermodynamic equation is converted to finite differences 
and made compatible with the Arakawa D-grid finite difference 
template developed for MODEL I (Achtemeier, et al . . 1986) . Fig. 1 
shows the template with the locations of the variables that appear 
in the thermodynamic equation. Note that the local tendency of 
temperature has been defined as the dependent and adjustable 
variable, E T . The finite difference version of the thermodynamic 
equation is, 


r 0 [e t + (m) xy(u) x °( T r ) y + (m) xy (v)y°( t v ) x 


+6 [ (Ip *y°-JL {ql) (T) ] ] +o 0 o 

c p 

-RjT t ) ^.4- - 2 --° 


C p F 


CgHF c E 


( 16 ) 


where the various overbar averages are defined in Achtemeier, et 
al. . (1986). 

3. MODEL II. Variational Equations 

The variational analysis melds data from various measurement 
systems at the second stage of a two-stage objective analysis. All 
data are gridded independently in the first stage and are combined 
in the second stage. The gridded observations to be modified are 
meshed with the dynamic constraints through Sasaki's (1970) 
variational formulation which requires the minimization of the 
integrand of an adjustment functional. Now it is not necessary to 
reproduce the full derivation of MODEL I plus the thermodynamic 
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equation in order to get MODEL II. Each term of the equations is 
a linear combination with the other terms. Therefore, all that is 
required is to perform the variational operations upon the 
thermodynamic equation and add the resulting terms to the 
appropriate adjustment equations of MODEL I. Let, 

I-2\ 5 m 5 +% 6 (E t -Et) 2 ( 17 ) 

where n a is the precision modulus weight for the temperature 
tendency and m 5 is equation (16) . Performing the variations upon 
each of the dependent variables that appear in the thermodynamic 
equation yields the following terms to be added to the respective 
variational equations. 


6E T -n e ( E T -Ef ) +i? 0 A 5-0 

( 18 ) 

8 u-R o m y a 5 (T x ) y)*° 

( 19 ) 

8 v-Rjm) x (\ 5 (T y ) x ) yo 

( 20 ) 

8 6-R 0 X s [ {T a ) xyo -— (^)^(T) +X 5 o 0 

( 21 ) 


p 


81*— R„[® '(i) X -R c [ (in) >*(v) '(Ip*’) y 

- R „ [ „ -R * [ ( 517 ) v <b* <^> 

C P 


( 22 ) 
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Table l summarizes the modifications of the existing MODEL I 
equations that are required to implement MODEL II. The first 
column labeled "variable referenced" locates the variable in the 
grid templet shown in Fig. l to which the new terms are referenced. 
For example, the new terms to be added to the existing function F 1 
(first line in Table 1) are calculated for the location of u in 
Fig. 1. Also included are two new equations, the latter being the 
thermodynamic equation. This brings to 13 the number of linear and 
nonlinear equations to be solved. 

4. MODEL II: Evaluation 

The purpose of this section is to demonstrate whether MODEL II 
performs as predicted by theory. In our evaluation of the 
variational assimilation models, we have used three criteria which 
have found use in the verification of diagnostic analyses 
(Krishnamurti, 1968; Achtemeier, 1975; Otto-Bliesner, et al. . 
1977) . These criteria are measures of, first, the extent to which 
the assimilated fields satisfy the dynamical constraints, second, 
the extent to which the assimilated fields depart from the 
observations, and third, the extent to which the assimilated fields 
are realistic as determined by pattern recognition. The last 
criterion requires that the signs, magnitudes, and patterns of the 
hypersensitive vertical velocity and local tendencies of the 
horizontal velocity components be physically consistent with 
respect to the larger scale weather systems. 
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The strong constraint formalism requires that the dynamical 
constraints; the nonlinear horizontal momentum equations, the 
hydrostatic equation, an integrated form of the continuity 
equation, and the thermodynamic equation be satisfied exactly (to 
within truncation) . Therefore, it is appropriate that the first 
evaluation of the variational model determine whether indeed the 
adjusted fields of meteorological variables are solutions of these 
physical equations. 

In solving the Euler-Lagrange equations, we substituted 
observed or previously adjusted variables into the nonlinear terms 
and other terms that are products with the Rossby number or are 
higher order terms and treated these terms as forcing functions. 
This approach made the linearized equations easier to solve but 
several cycles with the forcing terms updated with newly adjusted 
variables were required for the method to converge to a solution. 

The technique for determining whether the method converges to 
a solution is as follows. First, we note that any variable is 
found from the algebraic sums of all other terms of an equation. 
Thus the residual obtained by substituting variables back into the 
equation will be identically zero - the equation is satisfied 
exactly. This does not mean that the variational method has 
converged. Entirely different values for all of the variables may 
be found at the next cycle. Therefore, the adjusted variables are 
averaged over two successive cycles. Then the averaged variables 
are reintroduced into the dynamic constraints. Residuals are 
computed as remainders of algebraic sums of the terms of each 
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constraint. The root-mean-squares (RMS) of these differences 
(Glahn and Lowry, 1972) vanish when variables at two successive 
cycles are unchanged. When this happens, the constraints are 
satisfied and the method has been judged to converge to a solution. 
A convenient measure of how rapidly the method is converging to a 
solution is the percent, reduction of the initial unadjustment given 
by, 

Ar (%) -100 (1- - °~ rt ) (23) 
r ° 

The performance of MODEL II is assessed through the percentage 
reductions in the RMS differences from the initial unadjustments 
through the first four cycles of the solution sequence. The 
calculations are done for the eight adjustable levels in the model. 
Table 2 shows the percentages for the two nonlinear horizontal 
momentum equations. These results compare favorably with the MODEL 
I percentage residual reductions. The initial unadjustments are 
approximately halved at each cycle to about 90 percent after four 
cycles. 

The percentage reductions of the initial unadjustment for the 
integrated continuity and hydrostatic equations are shown in Table 
3. The RMS differences for the integrated continuity equation are 
reduced by from 96 to 99 percent at the second cycle and improve 
slowly to near 100 percent by ;the fourth cycle. These improvements 
are, of course, dependent upon the magnitudes of the initial 
unadjustment. We set the initial vertical velocity to zero. Then 
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the initial unadjustment is equal to the divergence integrated 
upward. The MODEL I cyclical solution order subjects the adjusted 
velocity components to a second adjustment to satisfy the 
integrated continuity equation. In this case, the averages of the 
adjusted velocity components are just averages of two solutions of 
the integrated continuity equation. Therefore the unadjustment 
should approach zero by the second cycle. 

The initial unadjustments for the hydrostatic equation at 
levels 4 through 8 are halved at each cycle and the percentage 
reduction increases to near 94 percent by the fourth cycle. 
Convergence is much slower at levels 1 and 2. There is a 65 percent 
reduction in the initial unadjustment at the second cycle at level 
2. There is no change during the third cycle and a slight increase 
in the initial unadjustment is observed at cycle 4. Given that the 
only difference between the adjustments presented here and the 
adjustments presented for MODEL I is the introduction of the fifth 
constraint, we are led to suspect that the degradation is directly 
related to the thermodynamic equation. 

Table 4 gives the percentage reductions of the initial 
unadjustment for the thermodynamic equation. Negative percentages 
occur where the RMS differences exceed the initial unadjustment. 
Table 4 shows that the initial unadjustment was reduced by nearly 
90 percent by the fourth cycle at levels 2 and 9. At the remaining 
levels, first cycle reductions of from 48 to 63 percent were 
followed by increases in the RMS differences that by the fourth 
cycle exceeded the initial unadjustment at levels 6 and 7. 
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Further analysis of the behavior of the convergence of MODEL 
II has revealed the following: 

1. The breakdown in the assimilation was almost exclusively in 
temperature. The initial unadjustments in the horizontal 
momentum equations and the continuity equation were reduced as 
was done with MODEL I. Only the first two levels in the 
hydrostatic equation showed any response to the temperature 
unadjustment and this was somewhat unexpected given that the 
most severe departures from convergence in the thermodynamic 
equation occurred at higher levels. 

2. The patterns of winds and heights generated by MODEL II (not 
shown) were unchanged from the winds and heights generated by 
MODEL I. The pattern analysis was an additional confirmation 
that the breakdown in convergence in MODEL II was largely 
confined to the thermodynamic equation. 

3. The initial unadjustment in the thermodynamic equation was 
found to be approximately an order of magnitude larger than 
the initial unadjustments for the other dynamic constraints 
and was approximately two orders of magnitude larger in the 
stratosphere. Although this is not the cause for the breakdown 
in convergence, it does show that a gross imbalance existed in 
the initial gridded fields of meteorological variables when 
those variables were substituted into the thermodynamic 
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equation. 

4. Analysis of the patterns of the residuals remaining after the 
fourth pass found that they were almost identical to the 
patterns of vertical velocity. 

Our analysis of the large RMS differences in the thermodynamic 
equation remaining after four cycles reveals the following 
concerning how the initial and adjusted vertical velocity adversely 
impacted upon the analyses. First, the initial vertical velocity 
was calculated kinematically and subjected to the variational 
adjustment by O'Brien (1970) . This method can transfer error from 
the lower levels into the upper levels of the troposphere and 
generate large and noisy vertical velocity patterns there. 
Furthermore, there is no consideration given for the change in 
static stability between the troposphere with its relatively large 
vertical velocities and the stratosphere with its relatively small 
vertical velocities. The kinematic vertical velocities were 
unrealistically large in the stratosphere and, when coupled with 
the large static stability, produced large and uncompensated terms 
in the thermodynamic equation. Therefore, the magnitudes of the 
initial unadjustments were approximately two orders of magnitude 
larger than were the initial unadjustments for the other dynamical 
constraints. 

Second, further theoretical analysis has revealed that the 
adjustment for the divergent part of the wind is the "weak link" in 
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this variational assimilation model. First order terms that contain 
the divergence adjustment cancel out in the cyclical solution 
formulations. The divergence adjustment must then be carried in 
second order terms and through other variables. Our solution for 
this problem has been to require the adjusted horizontal velocity 
components to satisfy the continuity equation constraint after each 
cycle, a variational model within a variational model, then allow 
for the second order terms and the readjusted velocity components 
to "nudge" the solution toward the desired dynamic balance. The 
result was that the RMS differences grew after the first cycle when 
the vertical velocity was released to converge slowly toward 
another equilibrium. 

5. Coupling the Vertical Velocity in MODEL I. 

In this section, we propose solutions for the vertical 
velocity related problems of very large initial unadjustments for 
the thermodynamic equation and the buildup of RMS differences in 
MODEL II. 

The solution for the problem of very large initial 
unadjustments in the thermodynamic equation is the implementation 
of a blended vertical velocity algorithm such as the variational 
method presented by Chance (1986) . This method, developed as part 
of this variational assimilation project but not included in the 
version of MODEL II evaluated as part of this study, blends the 
divergence of the horizontal wind with the vertical velocity 
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calculated from the adiabatic method. The relative weighting given 
the horizontal and the vertical velocity is a function of the 
stability, relative humidity, and satellite observed cloud cover. 
The divergence of the horizontal wind receives the greatest weight 
when the conditions of low stability, near saturation, or dense 
cloud cover at levels with near saturation prevail. The adiabatic 
vertical velocity receives greatest weight at locations where 
stability is high. Division by large stability reduces the 
magnitude of the vertical velocity in the stratosphere and forces 
the vertical velocity to near zero at the tropopause rather than at 
the arbitrarily defined top of the model domain. 

The formula for the modified vertical velocity is 


o M (k) 


d M (k-1) +D(k) Ao+ad T (k) 
l^a 


( 24 ) 


The modified vertical velocity at level k is the weighted sum of 
the modified vertical velocity at level k-1 plus the incremental 
vertical velocity obtained through the continuity equation and the 
vertical velocity obtained by the adiabatic method. The weight, a, 
carries the theoretical relative accuracies of the two methods for 
calculating vertical velocity as obtained through standard errors 
of observation for the observed variables. The weight also carries 
the relative 1 importance of the vertical velocities as determined by 
meteorological considerations. For example, the adiabatic 
vertical velocities are assigned the greatest weight in the 
stratosphere because the adiabatic method carries information 
regarding static stability, 
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’-‘-k-iry {T *-r T * ) 


(25) 


However, in the lowest layers of the analysis domain, a=0 to 
account for the near adiabatic conditions within the planetary 
boundary layer . 

Preliminary studies with the blended vertical velocity show 
that large magnitude centers of either sign developed by the 
kinematic method in the upper troposphere and lower stratosphere 
are reduced or eliminated. Therefore the large initial 
unadjustments that exist because of the use of the kinematic 
vertical velocities will be reduced or eliminated also. 

The solution for the problem of buildup of RMS differences in 
MODEL II is to reformulate the MODEL I variational equations so 
that the solution sequence will better couple the vertical velocity 
with the dynamic adjustment. Achtemeier, et al. (1986) have shown 
that the derivations in MODEL I required to reduce the number of 
dependent variables and equations to a single diagnostic equation 
in geopotential cancel out the zero order divergence adjustment 
terms. The adjustment of the divergent part of the wind is 
therefore forced into higher-order nonlinear terms which do not 
sufficiently impact upon the final adjustment to bring about 
compatibility with the continuity equation. The continuity 
equation was satisfied through the second variational step which 
forced an adjustment of the adjusted velocity components. The 
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problem was that the two variational steps could not be connected 
in a way that allowed adjustments required for satisfaction of the 
thermodynamic equation to feed back to the continuity equation. 

This analysis of MODEL II reveals that the second variational 
step must be eliminated and the coupling of the vertical velocity 
with the remainder of the adjusted variables must be part of a 
single variational model. It was found that the divergent part of 
the wind obtained from the first step adjustment. is a function of 
the nonlinear terms of the horizontal momentum equations. If F 5 

represents the nonlinear terms of the u-component equation and F 6 
represents the nonlinear terms of the v-component equation, then 
the horizontal momentum equations can be expressed as 

m i~- v+ -fa +F 5° 0 (26) 

m * mU+ -j£ +F *-° (27 > 

Forming the divergence from (23) and (24) and integrating through 
the depth of the analysis domain gives 

f(u x +v y ) do--f (F €x -F 5y ) do-0 (28) 

Equation (25) is an integral of the vorticity equation. The 
constraint upon the divergent part of the wind, and hence the 
vertical velocity, that must be satisfied in order for all MODEL I 
dynamic constraints to be satisfied is as follows. A particular 
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solution of the vorticity equation must integrate to zero at the 
top of the model domain - the particular solution being that the 
divergent component of the same adjusted wind field must also 
satisfy the integrated continuity equation. 

The incorporation of the integrated vorticity equation into 
the variational formalisms is the subject of MODEL IIB derived in 
Chapter III. 
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FIGURE CAPTIONS 

Fig. 1. The grid template for the variational assimilation model. 
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Table 2. Percent reduction of the initial unadjustment 
in the horizontal momentum equations after 4 cycles. 


Cycle Level 


No. 

2 

3 

4 

5 

6 

7 

8 

9 




u-component 




0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

54 

54 

52 

51 

50 

50 

51 

51 

2 

81 

78 

77 

75 

74 

75 

76 

76 

3 

92 

89 

87 

86 

86 

87 

87 

87 

4 

94 

93 

90 

89 

91 

91 

90 

90 




v-component 




0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

54 

53 

52 

53 

51 

51 

50 

50 

2 

78 

80 

77 

80 

- 77 

76 

76 

73 

3 

88 

89 

87 

90 

88 

88 

87 

84 

4 

93 

92 

91 

92 

91 

91 

91 

88 


Table 3. Percent reduction of the initial unadjustment 
in the integrated continuity and hydrostatic equations 
after 4 cycles. 


Cycle Level 


NO. 

2 

3 

4 5 

6 

7 

8 

9 

0 

0 

Integrated Continuity 
0 0 0 0 0 

0 

0 

1 

— 

— 

- — 

— 

- 

— 

— 

2 

97 

98 

98 99 

99 

99 

99 

99 

3 

96 

98 

98 99 

99 

99 

99 

99 

4 

96 

98 

99 99 

99 

99 

99 

99 

0 

0 

0 

Hydrostatic 
0 0 0 

0 

0 

0 

1 

51 

50 

50 50 

50 

50 

50 

50 

2 

73 

65 

75 75 

75 

75 

75 

75 

3 

83 

65 

88 88 

88 

88 

88 

88 

4 

86 

62 

94 94 

94 

94 

94 

94 
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Table 4. Percent reduction of the initial unadjustment 
in the thermodynamic equation after 4 cycles. 


Cycle 

No. 

2 

3 

4 

Level 

5 

6 

7 

8 

9 



Thermodynamic 

Equation 



0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

54 

60 

62 

63 

61 

63 

63 

48 

2 

81 

80 

74 

55 

24 

39 

76 

72 

3 

89 

73 

61 

32 - 

■12 

9 

62 

83 

4 

88 

65 

50 

14 - 

■38 

-12 

49 

89 
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Table 1. Modifications to variational equations in 
MODEL 1 to obtain MODEL 2. 


Variable Existing 

Referenced Function New Terms to be Added 


u 


Rjm) y(k 5 op y) xa 


R 0 (m) *u 5 (T y ) x )y° 


^o +R o^ [ ( t 0 ) «y*- -5- (qp ^(t) *y\ 

c p 


Eg 34 p 39 F 8 — {[ (m) x (u) *(Xp y0 ] x + [ (m) y (v) y (Xp 

Achtem. etal + [ (oX^)^ 0 ] 0 + — [ {o\ s ) xy q 2 + (w g l 5 ) } 

c p 

1986 


Eg 47 p 41 F 8 /y 
Achtem. etal 
1986 




New Equation X 5 (E T -Ef) 


d New Fgua fci on F r - -{ [ (m) *y ( u) x0 ( T x ) y + (m) xy (v) yo ( T y ) x ] 

+6 [ ( xya — — ( g ^) xy ( T ) xy ] +-^0 

c p ■*<> 

- (g 4 ) [ ^T r + (t) * y] ) 






'-2 - V 7 


preceding page blank not filmed 
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ABSTRACT 

A variational objective analysis technique that modifies 
observations of temperature, height, and wind on the cyclone scale 
to satisfy the five "primitive" model forecast equations is 
presented. This analysis method overcomes all of the problems that 
hindered previous versions - problems such as over-determination, 
time consistency, solution method, and constraint decoupling. A 
preliminary evaluation of the method shows that it converges 
rapidly, the divergent part of the wind is strongly coupled in the 
solution, fields of height and temperature are well-preserved, and 
derivative quantities such as vorticity and divergence are 
improved. Problem areas are systematic increases in the horizontal 
velocity components, and large magnitudes of the local tendencies 
of the horizontal velocity components. The preliminary evaluation 
makes note of these problems but detailed evaluations required to 
determine the origin of these problems await future research. 

1 . Introduction 

This study was designed to determine the feasibility of a 
constrained objective analysis based upon the variational 
methodology of Sasaki (1958, 1970) . The method uses as dynamic 
constraints the five primitive equations for a dry, adiabatic, and 
non-viscous atmosphere: the two nonlinear horizontal momentum 
equations, the continuity equation, the hydrostatic equation and 
the thermodynamic equation. The method is diagnostic, however 
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given the similarities between the dynamic constraints and the 
hydrodynamical equations of numerical prediction models, there 
exists a potential for extension of the technique to the derivation 
of initial fields for numerical models. 

The potential of the variational methods for multivariate 
objective analyses has been explored with many dynamic constraints. 
Some of the studies and the constraints used are: the geostrophic 
approximation (Sasaki, 1958), the continuity equation (0,Brien, 
1970; Dickerson, 1978; Sherman, 1978; Ray et al., 1978), divergence 
and vorticity (Schaefer and Doswell, 1979) , the balance equation 
(Stephens, 1970) , the two horizontal momentum equations (Lewis and 
Grason, 1972; Bloom, 1983), the two horizontal momentum and 
hydrostatic equation (Lewis, 1972), and the two horizontal 
momentum, thermodynamic, and hydrostatic equations (Achtemeier, 
1975) . 

Past attempts to develop a multivariate objective analysis 
based upon Sasaki's variational method with the five "primitive" 
equations as dynamical constraints have encountered several 
fundamental problems. Courant (1936) showed that the number of 
subsidiary conditions (dynamic constraints) must be at least one 
less than the number of adjustable dependent variables else the 
problem is overdetermined and a solution is not guaranteed. The 
over-specification problem must be solved as the five primitive 
equations form a closed set with five dependent variables. 

The Euler-Lagrange operations yield local tendencies of the 
Lagrange multipliers if the local tendencies of the temperature or 
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the horizontal velocity components are explicit in the dynamic 
constraints. Boundary conditions for these terms are unknown. The 
problem of time consistency in variational problems has been 
explored by Lewis (1980, 1982) and Lewis et al (1983). More 
recently, the time consistency problem has been found more 
tractable through use of the adjoint method (Lewis and Derber, 
1985; Talagrand and Courtier, 1987). 

Achtemeier (1975) found that the Euler-Lagrange equations 
decoupled the divergent part of the wind from the remainder of the 
adjustment with the result that the continuity equation was not 
satisfied. Attempts to constrain the local tendencies of velocity 
and temperature to require exact solution of the continuity 
equation did not solve the coupling problem (Achtemeier, 1979) . 

The methodology to circumvent the above problems and the 
theoretical development of a primitive equation variational 
objective analysis is presented in the next section (mathematical 
details are presented in Appendices A, B,and C.) The method is 
evaluated in Section 3 . 

2. Theoretical Development 

The objective analysis is designed for a terrain-following 
coordinate surface. We used a nonlinear vertical coordinate 
created from two functions that are piecewise continuous through 
the second derivatives. In this coordinate system, all coordinate 
surfaces above a reference pressure level are pressure surfaces. 
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The dynamical equations appear in their simplest form in pressure 
coordinates. Furthermore, hydrostatic truncation errors are 
confined to coordinate surfaces below the reference pressure level. 
The problems of reducing hydrostatic truncation error along 
terrain-following coordinate surfaces has been the subject of 
considerable investigation (Kurihara, 1968; Gary 1973; Sundqvist, 
1975, 1976; Janjic, 1977, 1989, and Achtemeier, 1990). The 
vertical coordinate is described in Appendix A.l. 

Subjecting the pressure gradient terms of the horizontal 
momentum equations written in terrain-following coordinates to the 
variational operations separates the two pressure gradient terms 
and combines the large, now uncompensated terms with terms from the 
other equations. These uncompensated terrain terms can dominate 
the adjustment. A test found that these terms generated large 
error that caused the variational method to diverge. 

The pressure gradient problem was solved by 
nondimensionalizing the dynamic constraints (Charney, 1948; 
Haltiner, 1971) and partitioning the hydrostatic terms to isolate 
the terrain part so that the variational adjustment could be 
performed on the meteorological partition. Appendix A. 2 presents 
details of this procedure. 

As regards the time consistency problem, Fjortoft (1952) found 
that the local change in the winds could be approximated by the 
translation of a weather system along an advective or steering 
current, usually a smoothed middle tropospheric wind. Therefore, 
the local tendencies of the velocity components were partitioned 
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into advective components, represented by the steady part of a 
weather system moving within a steering current, and developmental 
components, represented by the development of a weather system. 
Appendix A. 3 describes the partition. The developmental components 
of u and v were defined as dependent variables to be subjected to 
the variational adjustment. 

Appendix A. 4 gives the five dynamic constraints as modified. 
Abridged forms of these equations are as follows: 


M^—V+ 4> X +DTU+HAU+ VAU+ EXT (M x ) -0 ( 1 ) 

M z - U+< J) y +DTV+HAV+ VAV+EXT (M z ) - 0 ( 2 ) 

M 3 -Q(u x +v y ) Ao+(d-6 o )+EXT(M 3 )-0 (3) 

Mi-b'+yT+EXTMt) -0 (4) 

M 5 -LTT+HAT+VAT+WT+6o 0 +EXT(M 5 ) -0 (5) 


Conventional symbols are used. Abridged terms are defined as 
follows: 

DTU(V) = developmental component of local tendency of u or v. 
LTT = local tendency of T. 

HAU (V or T) = horizontal advection of u (v or T) relative to 
a moving weather system. 

VAU (V or T) = vertical advection of u (v or T) . 

WT = product of vertical velocity with perturbations of 
stability. 
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EXT (M i ) = extra terms that arise from any of the following 

sources: 

a) Lambert conformal map image projection, 

b) conversion into the nonlinear vertical coordinate, 

c) expansion of the Coriolis and/or map scale factors. 

Q = a normalized pressure thickness weight that arises from 
(b) above. For pressure levels above 700 mb, Q = 1. 

The fourth term on the right hand side of (5) is the product of the 
layer average static stability with the vertical velocity. 

These equations have been nondimensionalized and terms 
expressed in powers of the Rossby number. All terms identified by 
three letters (eg. , LTU or EXT) are higher order terms - either 
multiplied by the Rossby number or of order 0.1 or terms that 
involve unadjusted (observed) variables. 

Dependent variables are u, v, $, a, T, e u , and € v . The latter 
two variables are the developmental components of the local 
tendencies of u and v. This formulation leaves five constraints 
and seven variables to be adjusted. 

Following Achtemeier (1975) , a variational objective analysis 
was developed for adjustments of the seven dependent variables 
subject to exact satisfaction of the dynamic constraints (l)-(5). 
As expected, the addition of the two new dependent variables (the 
developmental components of u and v) was sufficient to overcome the 
over-specification problem. As regards the time consistency 
problem, recomposition of the local tendencies of u and v from the 
advective and developmental components yielded tendencies that 
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compared favorably with observed 3-h changes in u and v. However, 
the decoupling problem remained. Attempts to readjust the 
divergent part of the wind by requiring the adjusted horizontal 
velocity components to satisfy the continuity equation through a 
"variational adjustment within a variational adjustment" were 
unsuccessful in satisfying all five constraints. 

An analysis of the growth of the divergent part of the 
adjusted wind was performed to determine how the variational 
solution decoupled from the continuity equation. It was found that 
the divergent part of the wind is determined by adjustments through 
the higher order terms (HOT) of (1) and (2) . The divergent 
components can be made to satisfy the continuity constraint if 
these higher order terms are made to satisfy a particular solution 
of the vorticity theorem. Define 

F 5 -HOT(M x ) (6) 


F 6 -HOT(M 2 ) 


(7) 


so that, 


M x --v+$ X +F s - 0 


( 8 ) 


M 2 -U+4> y +F 6 -0 


( 9 ) 
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Forming the divergence, 

u x +v y + F 6x -F 5y -0 ( 10 ) 

The function that must integrate to zero is 

J (u x +v y ) da-J (F 6x -F 5y ) do-0 ( 11 ) 

for the vertical velocity to vanish at the top at the top of the 
domain. Therefore, (11) is a particular solution of the integrated 
vorticity theorem, the particular solution also requiring that the 
horizontal divergence integrate to zero, a requirement for 
satisfaction of the continuity equation. 

It is necessary to build (11) into the dynamic constraints if 
the decoupling problem is to be eliminated from the variational 
objective analysis. Define F s and F 6 as dependent variables and 
revise the dynamic constraints as follows: 

M 1 —F^DTU+HAU*VAU*EXT(M 1 )~0 ( 12 ) 

M 2 =F 6 +DTV+HAV+VAV+EXT(M 2 ) "0 ( 13 ) 

M 3~0(F 6x -F 5y ) Ao+ (6-6 a ) +EXT(M 3 ) -0 ( 14 ) 


M 4 -<t> a +YT+EXT(M 4 ) -0 


( 15 ) 
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Mt-LTT+HAT+VAT+WT+dOa+EXTiMs) =0 ( 16 ) 

M 6 -V+<b x +F 5 -0 ( 17 ) 

Mj-U+4> y +F 6 - 0 ( 18 ) 

The variational objective analysis is developed from these seven 
constraints. The nine adjustable variables include the original 
seven plus F 5 and F 6 . 

The dynamical constraints are written on centered differences 
on an Arakawa D-grid (Mesinger and Arakawa, 1976) . The finite 
difference operators and finite averaging operators are defined 

i 

following Anthes and Warner (1978) . The conversion of the 
constraints from differential form into finite differences is given 
in Appendix B. 

The gridded fields of meteorological data to be modified are 
meshed with the dynamical equations through Sasaki's (1970) 
variational operations. To simplify the derivations, the 
frictional terms in the horizontal momentum equations and the 
diabatic heating term in the thermodynamic equation were set to 
zero. 

The finite difference analog of the adjustment functional is, 

F- 

i 3 


(19) 
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The integrand, I. . is 

'iJ 


(u-u°) 2 +n 1 ( v-v°) 2 +ti 2 (d-d °) 2 
+n 3 (<J>-4> 0 ) 2 +ti 4 (T-T°) 2 +tc 5 (4> x -4>£) 2 
+ 1X 5 (4> y -4>y) 2 + * 6 2 + *7 (e u "0 2 (20) 

+n 7 (e v -e°) 2 +ti 8 (e T -e?) 2 +tc 9 (F 5 -F s °) 2 
+ 7i 9 (F 6 -F 6 0 ) 2 + 2j;A i M 1 


The weights, 7T, , are Gauss' precision moduli (Whittaker and 
Robinson, 1926) . The gridded initial variables (u°, v°, a 0 , §°, T°, 
e u °, v °, F 5 °, F 6 °) enter in a least squares formulation and receive n. 
according to their relative accuracies. The strong constraints to 
be satisfied exactly are introduced through the Lagrangian 
multipliers, A i . 

Objectively modified meteorological variables are determined 
by requiring the first variation on F to vanish. A necessary 
condition for the existence of a stationary set is that the 
functions are determined from the domain of admissible functions as 
solutions of the Euler-Lagrange equations. The variation is to be 
carried out at every point (r,s) within the grid. Thus, upon 
setting the weights a, = bj = 1 and differentiating the integrand 
(20) with respect to the arbitrary variable a r s , the Euler-Lagrange 
operator in finite differences is 

The Euler-Lagrange equations resulting from the operations 
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3l i.j _ T 3g i.j 

3a a *'J 3a 

^r, s r, s 




( 21 ) 


specified by (21) are given in Appendix C [(C.7)-(C.16)]. 
Including the seven dynamic constraints, these complete a closed 
set of 17 of linear and nonlinear partial differential and 
algebraic equations. Solutions are difficult to obtain by 
conventional methods. Achtemeier (1975) proposed a cyclical 
solution method that moves higher order terms and terms involving 
unadjusted (observed) variables into forcing functions. These 
forcing functions may be expressed with observed variables at the 
first cycle and with previously adjusted variables at higher 
cycles. Therefore the forcing functions are known at each cycle. 
This method of solution is valid for the latitudes and scales of 
motion for which the Rossby number is less than one. 


Use of the cyclical solution method yields the following set 
of linear Euler-Lagrange equations: 

M 3 m -Qs( F tx~ F s y) Ao+(d-d 0 ) +q- 5 F 7 Ao=0 (22) 

M 4 -$ g +yT+$-0 (23) 


M 6 — v+$ x +F 5 -0 


(24) 
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Mj-u+$ y +F 6 - 0 

+ n e 4> 0<f + n 5 x 4> * + ic 5y <j> y y + tc 6 0 <t>° 
-n 3 <J>+A 5x +A, 6y +A 4a +F 4 “0 

(25) 

(26) 

WjU+Xg+Fi-0 

(27) 

71iV-1 5 + F 2 = ° 

(28) 

tc 2 o+A 3 +F 3 -0 

(29) 

n 4 r+yX 4 +F 8 “0 

(30) 

i* 7 (e u -e°) +i? o A 1 =0 

(31) 

i* 7 (e v -e°) +i? o A 2 -0 

(32) 

Tig (e r -e?) +f? o A 7 -0 

(33) 

n 9 (F 5 -F 5 °) -A 1 +A 5 -Ao (g 5 If) y -0 

(34) 

1*9 ( F^~ F ( j°) -A 2 +A 6 +Ao (g 5 X°) *~0 

(35) 


As shown in Appendix C, variables may be easily eliminated among 
the equations. There results three diagnostic equations in 
geopotential, vorticity, and divergence, 


(Wf+n 


x ) V 2 <J)+ ( -^ ) 4) 00 +n 5x <|) ^+7r s <j> ^ ( n 6o + ) <j>"-7t 3 <|) 
*10 Y Y 2 


• -G 4 “ 




7C 


10 


(®2y-®3x) 


( 36 ) 
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*ioC-*n V2 4>-G 2y ~G 3x (37) 

V 2 [ (&o) 2 q£n 2 D] -n 10 D-'V 2 G 1 +G Zx +G 3y (38) 

Details regarding symbol definition are found in Appendix C. 

The variational theory specifies natural boundary conditions 
that are consistent with the Euler-Lagrange equations. If it is 
assumed that there are no adjustments in the data along the 
boundaries, then the boundary conditions may be specified. In the 
latter case, the Lagrange multipliers, A.., are zero at the 
boundaries and the initial unadjusted values are used for the 
boundary conditions. 

Initially, the Euler-Lagrange equations were solved with 
specified boundary conditions. These boundary conditions forced 
high frequency waves into the solutions for the velocity components 
near the boundaries. Divergences calculated from these velocity 
components gave large erroneous vertical velocities. We therefore 
returned to the natural boundary conditions. 

The Euler-Lagrange operator for natural boundary conditions 

is, 



( 39 ) 
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Performing the operation specified by (39) produces a set of Euler- 
Lagrange boundary eguations in $, u, v, and D. Details of the 
derivations are given in Appendix D. The boundary conditions for 
§ are, 


( W + ( u°-F 5 ° ) -0 

5 Ttj + ttg ^ 5 ^l + ^9 




k x +tz 9 


) 4 > y - rc s 4 > y + J * 1 ” 9 ( u °+0 "0 


(40) 


The boundary conditions for u and v that are consistent with (40) 
are, 


(* 5 +*i 


tt 5 + Tl e, 




* 5 + * 9 


) V-* 5 (<J>£-F 5 °) — V °“0 

7l 9 

) u+n 5 (Qy+F 6 °) -n 1 715+719 u °-0 


(41) 


The derivation of (40) placed a constraint upon the boundary 
conditions for the divergence, namely, that the divergence must be 
specified along two rows or columns at the boundaries. 

Subject to the boundary conditions and specification of the 
precision moduli, (36) -(38) may be solved for the geopotential, 
vorticity and divergence. Coefficients for the second order 
partial derivative terms are always positive, the equations are 
elliptic, and thus solutions by standard methods are assured. Then 
u and v must be retrieved from the vorticity and the divergence. 

A number of investigators (Sangster, 1960; Hawkins and 
Rosenthal, 1965; Shukla and Saha, 1974; Schaefer and Doswell, 1979; 
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Lynch, 1988) have proposed methods for reconstruction of the 
velocity components from the vorticity and divergence (or 
streamfunction and velocity potential) . After investigating 
several of these methods, including those of Endlich (1967) and 
Bjilsma et al. (1986), it was determined that the Lynch method 
could be best adapted to the Arakawa D-grid with a minimum of error 
in reconstructing the velocity components. 

First, the field of divergence was modified by a small 
constant so that Gauss' theorem, 

| fDda-j V n ds ( 42 ) 


was satisfied. Then the u-component was reconstructed through 

(43) 

subject to mixed boundary conditions in u (obtained from (41)) 
along the y-boundaries and u y obtained along the x-boundaries from 

u y“ v x~Z (44) 

Then, beginning at the lower x-boundary with v from (41) , 

Vy‘D~U X 


(45) 
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was solved to find v uniquely. 

3. Case Study Description and Preprocessing of Data 

The data used to test the variational objective analysis 
consisted of rawinsonde temperature, height, and wind data at 
standard National Weather Service reporting sites shown in Fig. 1 
for a large part of the United States and parts of southern Canada 
on 12 GMT 10 April 1979 and 00 GMT 11 April 1979. This case was 
originally selected because microwave temperature soundings 
coexisted with special 3 -hr rawinsonde data over a large area of 
the central United States (small dashed-line box in Fig. 1) during 
a major cyclogenesis. The 3 -hr rawinsonde data were used as ground 
truth for the local tendencies of the velocity components and 
temperature diagnosed from the variational objective analysis. 

The data at 12 GMT 10 April 1979 described a weak, dissipating 
short wave moving northward over the Central Plains in advance of 
a more vigorous short wave. At 00 GMT 11 April, an intense jet 
streak moved northeastward over Oklahoma and Texas and triggered a 
mesoscale convective system over northern Texas that produced a 
number of fatalities at Wichita Falls, Texas. 

The data were gridded from the observations by a modification 
of the Barnes (1964) objective technique that is designed to 
minimize analysis error at the boundaries of the field of data 
(Achtemeier, 1986) and to provide accurate derivatives within the 
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interior of the domain (Achtemeier, 1989) . The analyses were done 
for 10 levels from 1000 mb to 100 mb. The horizontal grid was a 
40x25 array with a 100 km grid spacing. Then thermodynamic data 
were converted to the nonlinear vertical coordinate through a 
hydrostatically consistent interpolation downward from the 
reference pressure level of 700 mb to the terrain-following 
coordinate surfaces. In addition, a smoothed version of the 600 mb 
wind velocity components was obtained through a single pass of the 
objective interpolation designed to reproduce the long wavelength 
features inherent in the data. The smoothed wind field served as 
the advective wind in the calculation of the advective part of the 
local tendencies of the velocity components. 

The above analyses produced gridded fields of temperature, 
height, and u and v wind components. The initial fields of 
vertical velocity, developmental components of the local tendencies 
and F 5 and F 6 must be estimated from these data. Letting 


e 




dx dy 


(46) 


the adiabatic vertical velocity can be found by solving (B.10) for 
a. Then an adjusted vertical velocity can be found by a 
variational formulation using the continuity equation (Chance, 
1986) that is similar to the O'Brien (1970) method with the 
exception that compatibility between the divergence and the 
vertical velocity is forced at each level. The relative weight 
accorded to the adiabatic vertical velocity is directly 
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proportional to the static stability. Thus the adiabatic vertical 
velocity receives the greater weight in areas of higher stability 
such as the stratosphere. This procedure keeps large erroneous 
vertical velocities generated by divergence error from being 
transferred from the troposphere into the stratosphere where, in 
product with the static stability terms of (B.10), would produce 
large errors in the adjusted time derivatives of temperature. 

We have no way of estimating the developmental components of 
the velocity component tendencies from data collected at a single 
time. Therefore, these fields were set to zero. An alternative, 
if available, would be local tendencies from a numerical model. 

The forcing function variables, F 5 and F 6 are estimated by 
substituting the initial variables into (B.4) and (B.6). Then F 5 
and F 6 were adjusted to satisfy (11) with the exception that the 
integral of the divergence was replaced by the adjusted vertical 
velocity. 

The resulting fields (and selected derivative fields) of T, §, 
u, v, a, e u , e v , F 5 , and F 6 were designated as unadjusted fields and 
entered into the variational objective analysis through the 
functional integrand, I, given by (20) . The unadjusted quantities 
were accorded precision modulus weights according to the formula. 


Tt 




Gj (x, y) 

2a\ 


(47) 
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where the a i is the root-mean-square (RMS) error of observation for 
the ith variable. G,- is in general a function of observation 
density but G=1 for this study. However, since observational 
errors are available only for u, v, #, and T, only n 1f 7r 3 , and 7r 4 
can be obtained from (47) . the a i for the remaining unadjusted 
quantities must be inferred from the known observational errors 
through the dynamic constraints or simplifications therefrom. 
These a. are given by, 
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Here S is the average separation between observation sites. 

In addition, n 9 = as terms such as the Taylor series 
expansion of the Coriolis parameter in product with the wind are 
considered equal in weight with the wind itself. 

Table 1 shows the standard errors of observation for the 
winds, heights, and temperatures and the RMS errors for the other 
adjustable meteorological variables. Estimates for the scaler wind 
speed as functions of elevation angle of the balloon (Fuelberg, 
1974) are given in the first two columns. The values for the 20 
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degree elevation angle compare favorably with the results from 
Hovermale's (1962) spectral decomposition of meteorological data. 
RMS values for heights and rawinsonde temperatures are from a 
composite of methods for estimating measurement error (Achtemeier, 
1972) . 

Table 2 gives the nondimensional precision modulus weights 
calculated from the various functional relationships of the RMS 
errors from Table 1. The more accurately measured (estimated) 
variables receive larger values. Largest weights are accorded the 
geopotential height followed by the winds and temperatures. The 
developmental components of the local velocity tendencies receive 
the smallest weights. 

Several modifications in the n. given in Table 2 were made 
before the April 10-11 data were subjected to the variational 
objective analysis. First, the precision modulus weights for 
levels 9 and 10 of the vertical velocity were assigned large values 
to require the adjusted vertical velocity to vanish at the top of 
the domain. Second, the weights for the geopotential were reduced 
by a factor of 10 because prior studies gave solutions that were 
forced too strongly toward the geopotential. It is possible that, 
as a boundary condition, the geopotential has a greater impact upon 
the the solution than suggested by the magnitude of its precision 
modulus weight. 
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4. Evaluation of the Variational Assimilation Model 

Three diagnostic criteria were used to evaluate the 
variational objective analysis. These criteria are, satisfaction 
of dynamical constraints, adjustment departures from observations, 
and pattern analysis. 

a) Satisfaction of Dynamical Constraints 

The method must converge regardless of how well the other 
criteria are satisfied. But some method must be developed that 
demonstrates that the analysis does converge. The Sasaki (1970) 
strong constraint formalism requires that the dynamical 
constraints; the nonlinear horizontal momentum equations, the 
hydrostatic equation, the continuity equation, and the 
thermodynamic equation be satisfied exactly (to within truncation) . 
Recall that the cyclical solution method for solving the Euler- 
Lagrange equations required the substitution of observed or 
previously adjusted variables into the forcing functions. As a 
measure of progress toward convergence, at the end of each cycle, 
the adjusted variables were averaged with their respective values 
at the previous adjustment, reintroduced into the dynamical 
constraints and residuals calculated. It follows that the 
residuals decrease as the differences between adjusted variables at 
two successive cycles decrease. The residuals vanish (the 
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variational objective analysis converges) if the adjusted variables 
at two successive cycles are the same. A convenient measure of how 
rapidly the method is converging to a solution is the percent 
reduction of the initial unadjustment given by, 

Ar (%) -100 (1- - °~ rT ) ( 49 ) 
r ° 

Fig. 2 shows how the reductions of the initial RMS differences 
for the two horizontal momentum equations varies for each pass 
through the cyclical solution sequence for the eight adjustable 
levels of the model. The residuals for the u-component momentum 
equation are approximately halved with each cycle through the 
fourth cycle. The solution stabilizes to near 99-100 percent 
reduction of the initial unadjustment except for a 97 percent 
reduction at the 9th level after eight cycles. The RMS differences 
for the v-component equation decrease at the first cycle and level 
off at the second cycle. These differences increase slightly at 
level 7. Then the residuals decrease monotonically through the 
eighth cycle with reductions of the initial unadjustment of from 
98-99 percent (96 percent at level 9) . 

There were two reasons why the analysis was done through eight 
cycles. First, the objective of obtaining near 100 percent 
reduction in the RMS differences was accomplished for most levels. 
Second, regardless of the care taken in formulating consistent 
boundary conditions, there remained deleterious boundary effects 
that were drawn into the interior of the domain one grid space for 
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each cycle. The outer three rows of grid points were deleted from 
the evaluation statistics (see large dashed rectangle in Fig. 1) . 
Therefore, the effects of the boundary conditions entered the 
evaluation area beginning at the fourth cycle. 

The reductions of the initial unadjustment for the integrated 
continuity equation are shown in the left panel of Fig. 3. The 
rate of percentage reductions drops off after a large decline at 
the first pass but still reductions by the eighth pass were mostly 
between 97-99 percent. The slower convergence at level 9 (92 
percent after 8 cycles) and also at level 9 for the horizontal 
momentum equations may have been the result of large adjustments of 
the divergent part of the wind required for mass consistency with 
small vertical velocities in the stratosphere. 

The initial unadjustments for the hydrostatic and 
thermodynamic equations (middle and right panels of Fig. 3) 
monotonically decreased by about one half at each cycle. The 
percentage reductions of the RMS differences were mostly near 100 
percent at all levels by the eighth cycle. 

The satisfaction of constraints test shows that convergence 
toward a solution was obtained for all levels and for all five 
dynamic constraints. Therefore, MODEL IIB represents a significant 
advancement over the MODEL II. 
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b) Adjustment Departures from Observations 

The transferral of the observations to the grid and the 
modification of the gridded data to satisfy the dynamical 
constraints is a two-step process. Information from the 
observations is not available to the second step. Therefore, there 
is an implicit assumption that the initial gridded fields correctly 
carry the phenomena described by the observations. This assumption 
is not strictly true and it is necessary to grid the data with 
sufficient accuracy so that analysis error does not dominate the 
first and second derivatives. We have modified the widely used 
Barnes (1964, 1973) method for gridding meteorological data to 
yield significant improvement in the accuracy of the gridded data 
and its derivatives (Achtemeier, 1986, 1989) . 

In the section under a) above, we showed that the variational 
objective analysis converges to a solution. Now we seek to find 
whether the variational method improved upon the unadjusted 
analysis by adjusting the fields to better fit the original 
observations . 

Consider an "accuracy index" given by the solid lines in Fig. 
4. We first calculated two sets of RMS differences, one between 
values from the unadjusted fields at observation locations and the 
observations and the second between the adjusted fields and the 
observations. Then we subtracted from these RMS differences the 
standard errors of observation for wind components, height, and 
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temperature listed in Table 1. This means that if the results are 
zero, the objective analysis has gridded the data to within the 
standard error of observation for the data. If the results are 
negative, then the objective analysis has produced a better fit to 
the observations. Positive values mean that the adjustments have, 
on the whole, departed farther from the observations than expected. 
In interpreting these results, it must be kept in mind that the 
mean winter standard observational error estimates taken from 
Hovermale's (1962) results do not exactly express the true 
observational error for this case. Thus, some small departure of 
either sign from given values should be expected. 

The accuracy index for the unadjusted and adjusted heights and 
temperatures (Figs. 4a and 4b) were within acceptable limits. The 
index for the adjusted heights was displaced toward the positive, 
an indication that adjustments away from the observations were 
necessary to bring the fields into constraint satisfaction. The 
unadjusted fields of the horizontal velocity components were also 
within acceptable limits (Figs. 4c and 4d) . However, above 800 mb, 
large positive values for the adjusted velocity components show 
that the variational analysis produced wind fields that were, 
significantly different from the observations. 

The dashed lines in Fig. 4 are the means of the differences 
between the unadjusted (adjusted) fields interpolated to the 
observation sites and the observations. Means near zero are 
expected unless systematic adjustment is required to > achieve 
solution of the variational equations. Means were near zero for 
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the heights and the temperatures, except for temperatures near the 
tropopause between 300 mb and 200 mb where systematic adjustments 
were expected. The means were also near zero for the unadjusted 
velocity components. However, large systematic adjustments were 
found for the variationally adjusted velocity components (Figs. 4c 
and 4d) . The u-components were increased on the average 6 m s' 1 
between 500 mb and 300 mb. The v-component systematic reduction 
was a linear function of pressure. The v-component was on the 
average decreased by approximately 8 m s' 1 between 300 mb and 200 
mb. 

There was no systematic modification of the height fields that 
could be called upon to explain the adjustments in the velocity 
fields. An error in the mathematical derivation of the dynamic 
constraints or in the programming is suspected in these cases. The 
pattern analysis should provide further insight into the origin of 
these large systematic adjustments. 



63 


c) Pattern Analysis: 00 GMT 11 April 1979 

Maps of heights, wind vectors, and temperatures were taken 
from selected levels within the domain of the variational objective 
analysis for 00 GMT, 11 April 1979, in order to interpret the 
statistical results presented in subsections a) and b) . 
Comparisons were made between patterns in the unadjusted initial 
fields and the adjusted fields. The analyses were done on the 
synoptic scale however, we note that a mesoscale convective system 
was located within the high frequency observation area over parts 
of Texas and Oklahoma. 

Heights at 60 m intervals and wind vectors at 300 km intervals 
are shown in Fig. 5 for 800 mb, 500 mb, and 300 mb. The convention 
for wind speed is: flag (25 m s' 1 ), barb (5 m s' 1 ), and short barb 
(2.5 m s' 1 ). At 800 mb, the circulation center has been displaced 
from its unadjusted location over northwestern Colorado to its 
adjusted location over eastern Colorado in better agreement with 
the center of lowest heights. Elsewhere, adjustments in both 
heights and winds at 800 mb were small (Fig. 6) . At 500 mb (Fig. 
5) , the unadjusted analysis placed a weak short wave trough 
oriented eastward into Kansas from the parent trough. No trough 
appears in the wind field over Kansas. Thus, winds blow from high 
to low heights over Texas and Oklahoma and from low to high heights 
over Nebraska. The adjusted winds have been turned to more 
westerly in better agreement with the heights over Texas and 
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Oklahoma however, east of the Great Plains, the adjusted winds turn 
to blow toward higher heights. The same pattern of adjustment is 
also evident at the 300 mb jet stream level. The unadjusted 
analysis fits the winds with the height field. The adjusted 
analysis increases the wind speeds and turns the winds more 
westerly to blow toward higher heights. 

The differences between the adjusted and unadjusted fields are 
shown in Fig. 6 for 500 mb and 300 mb. In general, the variational 
objective analysis decreased the heights over the northern states 
and increased the heights over the southern states. The 10 m 
adjustment over Oklahoma at 500 mb tended to lessen the sharpness 
of the short wave trough there. Elsewhere, heights were lowered 
15-20 m over Montana. 

Fig. 6 also shows that an average 5ms' 1 westerly component 
was added to the wind field at 500 mb and an average 10 m s' 1 
northwesterly component was added to the 300 mb wind field. This 
broad scale adjustment has no apparent relationship to either the 
height field adjustment or the synoptic weather pattern. 

Fig. 7 shows fields of unadjusted and adjusted mean layer 
temperatures for 750 mb, 450 mb, and 250 mb. The unadjusted 
patterns at all levels have been preserved by the variational 
objective analysis. Temperature adjustments were less than one 
degree at 750 mb and 450 mb. The variational analysis cooled the 
250 mb layer by an average of 2C. The unadjusted layer average 
temperature was too warm across the tropopause and the change was 
made to make the temperatures consistent with the heights. 
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The variational objective analysis modified height, 
temperature, and wind velocity for satisfaction of the dynamic 
constraints. We now assess how these adjustments have impacted 
upon derivative quantities such as vorticity, divergence, and 
vertical velocity that are derived from these basic fields. In 
addition, the local tendencies of the velocity components and 
temperature are determined from arithmetic sums of adjusted terms. 
Patterns of these sensitive variables must be physically realistic 
when compared with other data sets such as cloud fields, 
precipitation, and independent measurements of the variable. 

Patterns of relative vorticity for the unadjusted and adjusted 
wind fields are shown in Fig. 8 for 500 mb. The variational 
objective analysis shifted the vorticity gradient, identifying the 
area of positive vorticity advection and upward vertical velocity, 
from over the Texas panhandle to over Oklahoma and Kansas, 
locations coincident with the mesoscale convective system. 
Elsewhere, there were only small differences between the unadjusted 
and adjusted vorticities. 

A comparison of the 500 mb vertical velocity patterns (Fig. 9) 
shows that the variational objective analysis shifted the center of 
maximum vertical velocity eastward from the Texas panhandle to 
western Oklahoma in better agreement with the location of the 
mesoscale convective system located over central Oklahoma and north 
Texas. The variational analysis also weakened the subsidence area / 
over Nebraska by 2 cm s' 1 . The subsidence area over Louisiana and 
eastern Texas in the unadjusted vertical velocities was replaced by 
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2-4 cm s' 1 rising motion in the adjusted field. Deep convective 
precipitation was located within this area (see shaded area in Fig. 
9) • 

Once the variational objective analysis was completed, the 
developmental components of the local tendencies of velocity 
components and temperature were recombined with the advective 
components, redimensionalized, and expressed as 3-hr changes. 
These 3-hr "adjustment" tendencies were compared with tendencies 
calculated from 3-hr rawinsonde data collected over the central 
part of the United States as part of the NASA AVE/ SESAME project 
(see fine dashed grid in Fig. 1) . Then "unadjusted" 3-hr 
tendencies were calculated upon substitution of unadjusted 
variables into the dynamical constraints and solving for the 
tendency terms. Inherent in these comparisons is an assumption 
that the observed 3-hr tendencies are "ground truth". This 
assumption is not strictly valid for the following reasons. First, 
it is likely that some of the observations, either at 0000 GMT or 
at 0300 GMT, were influenced by the mesoscale phenomena within the 
analysis areas. Second, the unadjusted and adjusted 3-hr 
tendencies were calculated from 0000 GMT data and are therefore 
centered at 0000 GMT. These tendencies were compared with the 
ground truth tendencies that were calculated from observations 
taken at both 0000-0300 GMT and are therefore centered at 0130 GMT. 
And third, extrapolation of the local tendencies calculated from 
the unadjusted and adjusted data has validity only if the time 
scales of the passage of the weather systems are much greater than 
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three hours. 

Fig. 10 shows fields of the 3-hr u-component tendencies at 800 
mb and 500 mb. The observed tendencies show increases in u over 
Oklahoma and decreases in u over northern Missouri and Iowa. Both 
unadjusted and adjusted tendencies show similar features but they 
are shifted to the southwest by about 500 km. Note also that the 
unadjusted and adjusted tendencies have approximately the same 
pattern and the centers from the variational objective analysis 
tend to be slightly larger in magnitude. 

The v-component tendencies at 800 mb and 500 mb are shown in 
Fig 11. Unlike the u-component tendencies, the centers for 
unadjusted and adjusted tendencies are approximately collocated 
with the observed centers. The magnitudes of the positive center 
over Arkansas compare well at 800 mb however the adjusted field 
shows little correlation with the observed v-tendencies in the 
western half of the grid. At 500 mb, the centers were mostly 
collocated however the magnitudes for both the unadjusted and 
adjusted v-component tendencies were much greater than the observed 
3-hr magnitudes - the magnitudes of the adjusted v-component being 
the largest. 

At 300 mb, Fig. 12, both unadjusted velocity component 
tendencies departed considerably from the observed fields. The 
adjusted tendencies appeared to be no more correct. 

Table 3 gives correlation coefficients between the unadjusted 
(initial) and observed tendencies and between the adjusted 
(variational) and observed tendencies for the eight interior levels 
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of the analysis domain. Somewhat surprisingly, the adjusted 
correlations were higher than the unadjusted correlations for most 
levels below 500 mb. In calculating the correlation coefficients 
that appear in Table 4, we shifted the adjusted and unadjusted 
tendency fields to the northeast approximately 150 km to account 
for the 1.5 hr translation of the weather system. The correlations 
for the shifted tendencies were larger. The variational objective 
analysis gave improvement over the unadjusted u and v tendencies 
however, in general the correlations for the adjusted fields were 
in the range from 0.5-0. 8 below 500 mb and were still negative 
above 400 mb. Results for the temperature tendencies in both 
tables showed no clear indication of superiority of the adjusted 
temperatures over the unadjusted temperatures. 

5. Discussion 

Based upon our experience with developing a basic variational 
objective analysis technique (Achtemeier et al., 1986) we have 
derived a new variational objective analysis method that appears to 
solve all of the problems encountered with earlier versions. These 
problems included the problem of over-determination noted by 
Courant (1936), the problem of time consistency that arose upon 
applying the direct variational method to local tendencies of wind 
velocity components and temperature, the problem of solving a set 
of complicated nonlinear partial differential equations, and the 
problem of decoupling the divergence equation constraint from the 
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remaining dynamical constraints. This version of the objective 
analysis contains more equations and requires more complicated 
solution methods than were necessary for the 1986 version. 

The evaluation presented in this report is only preliminary in 
that it identifies problems with the method but does not determine 
whether the problems are endemic to the method and therefore 
degrade data assimilation or whether the problems arise because of 
correctable errors in the mathematical derivations or the 
programming. 

The satisfactory results of the evaluation are as follows. 

1) The method converges for all five dynamic constraints. The 
divergent part of the wind is strongly coupled in the 
solution. Convergence after only eight cycles ranged mostly 
between 98-100 percent of the initial unadjustment with the 
poorest convergence at the 9th level still at an acceptable 92 
percent. 

2) The method gave reasonable adjusted fields of heights and 
temperatures from the standpoint of pattern recognition. The 
major synoptic weather systems were retained from an accurate 
initial objective interpolation to the analysis grid. Smaller 
features such as short waves were also retained. The method 
did not introduce erroneous wavelengths into the adjusted 
fields. 

3) Sensitive derivative fields such as vorticity and vertical 
velocity were better located with respect to important 
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precipitation producing weather systems relative to the 
unadjusted fields. Gradients of positive vorticity advection 
were collocated with upward vertical velocity centers. 

The unsatisfactory results from the evaluation are as follows. 

1) The variational objective analysis systematically increased 
the zonal component of the wind in a way that caused 
significant departures from the original observations. These 
departures appeared to be a function of elevation and of 
latitude from the grid origin (the largest increases were 
found in the eastern part of the grid) . These departures 
systematically turned the winds east of the Great Plains to 
blow from low to high heights. 

2) Though at many levels, the patterns were similar, the 
variational objective analysis greatly overestimated the 
magnitudes of the local tendencies of the wind components and 
temperature. Correlations between verification 3 -hr 
tendencies and 3 -hr tendencies derived from adjusted data 
ranged from about 0.5 to 0.8 at levels below 500 mb. 
Correlations were mostly very small or negative at 200 mb and 
300 mb. 

The reasons for the unsatisfactory results await a more 
thorough analysis of the method. The systematic increases in the 
adjusted wind velocity are suggestive of an error embedded within 
the mathematical formulas or coding of the programs. We were able 
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to trace the vary large magnitudes of the tendencies to the 
advective components. These are relative simple formulations and 
it has yet to be determined why large advective changes in velocity 
were found in both the unadjusted and adjusted fields but were not 
observed . 

It could be argued that the large tendencies of the adjusted 
fields should have been expected given that a mesoscale convective 
system was within the analysis area during the period 0000-0300 
GMT. The variational objective analysis was rerun for 1200 GMT 10 
April 1979 data set to test this argument. This period was 
characterized by the same general synoptic scale long wave trough 
over the western United States. There were no significant 
precipitation systems active however. The results showed large 
magnitude centers of the local tendencies of u and v in both the 
unadjusted and adjusted fields. Therefore, the finding of large 
magnitude tendencies within the 0000 GMT 11 April variational 
analysis was not coincidental with severe weather. 

In conclusion, the variational objective analysis represents 
a mammoth effort in mathematical development and programming. One 
must question whether, if the problems encountered thus far are 
solved, the difficulty of the method would limit its use in routine 
analysis of meteorological data given that there are other 
nonvariational techniques for blending meteorological data that are 
being used with success. The answer to the question will in part 
be delayed until the methods currently in use have been fully 
applied and evaluated. 
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Appendix A: The Dynamic Constraints 


Following Shuman and Hovermale (1968) , the horizontal momentum 
equations and the continuity equation that form the basis of the 
numerical variational objective analysis/assimilation method are 
written below as they appear in an arbitrary vertical coordinate 
and cartesian on a conformal projection of the earth: 


du / du du \ • du ^ , dd> dp 3d> » ^ 

dt dx dy do dp dx dx 


(A.l) 


0(A.2) 

at ox dy do dp dy dy 


dt do dx dy do dx dy 


The hydrostatic equation is, 


da 3<t> + RT = q 
dp do p 


and the thermodynamic equation, 
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These equations must be subject to several transformations before 
they can be used in a successful variational method. These 
transformations are described in the following sections. 

A. 1 A Nonlinear Vertical Coordinate System 

The vertical coordinate is designed to concentrate horizontal 
variations with the lower coordinate surface to levels below a 
reference pressure level p*. The coordinate surfaces above p* are 
constant pressure surfaces. The transformation into a nonlinear 
vertical coordinate was done for the following reasons: 

(1) The dynamical equations appear in their simplest form on 
pressure surfaces. The complex, compensatory terms are 
confined to levels below p*. 

(2) Vertical interpolation of meteorological observations to 
coordinate surfaces is not required for pressure surfaces. 
Further, there is no need to interpolate from sigma 
coordinates back to pressure surfaces for purposes of 
interpretation of the variationally adjusted fields of data. 

(3) Hydrostatic truncation error and pressure gradient force 
errors are eliminated on the pressure levels above p*. The 
problems of reducing hydrostatic truncation error along 
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sloping coordinate surfaces are well known (Achtemeier, 1990) . 


Two curves that are piece wise continuous through the second 
derivatives make up the nonlinear vertical coordinate. The upper 
layer relates to pressure by a straight line. Boundary conditions 
are a = 0 at p = p u and a = a* at p = p*. This equation is, 


o-o* 


P~P U 
P*-P u 


(A. 6) 


Boundary conditions for the lower curve are a - 1.0 at p = p s and 


o-o* 


da o* 
dp~ (P,-P„) 


^o 

dp 2 
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(A. 7) 


at p = p*. The lower curve, a cubic polynomial, is, 


o-p (p-p*) 3 +o* — — , 
P*-P u 


(A. 8) 


where 


p.(l- 0 *£ g- Pu ) (p s -p*)~ 3 . (A. 9) 

P ~Pu 


Fig. A. 1 shows the distribution of coordinate surfaces below 
600 mb for the approximate range of surface pressures (800 to 1025 
mb) for the smoothed orography of the variational analysis. The 
reference pressure p* is 700 mb. These coordinate surfaces tend to 
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follow constant pressure surfaces at locations away from areas of 
high elevation. The compression of the coordinate surfaces over 
higher elevation is clearly evident. 

Other variables that are an outcome of the nonlinear vertical 
coordinate appear elsewhere in the transformation of the dynamic 
equations. These are: 


where , 



(A. 10) 


(p-p*) 2 +a, 

J s -3p (p 3 -p*) 2 +a . 

It is understood that if p - p* < 0, then p - p* a 0. 

Terms in the dynamic equations that must be transformed are as 


follows: 
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1 1 00 


Fig. A. 1 
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(1) The pressure gradient force terms of the horizontal momentum 
equations (A.l and A. 2) take the form, 


dp dx ox ox ox 


(A. 11) 

(2) The first term of the continuity equation transforms into 




(3) The hydrostatic equation transforms to, 


(A. 12) 


04 > T ainp _ Q 

do do 


(A. 13) 


(4) The fourth term of the thermodynamic equation (5) becomes, 


RTto RT / > .. \ 

TT” — (Q-3® + <34»,) 


(A. 14) 


A. 2 Reduction of Terrain Impacts upon Analysis 

Small hydrostatic residuals and related pressure gradient 
force errors that plague numerical models written in terrain- 
following coordinates have been well documented. Much larger 
errors can be generated upon subjecting the pressure gradient terms 
of the horizontal momentum equations to the variational operations. 
The variational operator separates the two pressure gradient terms 
and combines the large now uncompensated terms with terms from the 
other equations. The terrain terms, for which the 
nonmeteorological part may exceed 90 percent of the magnitude of 



79 


the term, can dominate the adjustment. A test found that these 
terms generated large error that caused the variational method to 
diverge. 

The above problem may be avoided if the hydrostatic terms are 
partitioned to isolate the terrain part so that the variation can 
be applied to only the meteorological "signal". Note that a 
partition not a transforation is done. There is no change in the 
vertical coordinate. 

The equations were nondimensionalized following the 
methodology of Charney (1948) and Haltiner (1971). The resulting 
nondimensional variables contain the "whole" signal. The 
geopotential height and temperature are partitioned into terrain, 
reference, meteorological, and residual categories according to, 

( ) „-( ) r +( ) *+< ) °+[< )-( ) °] (A. 15) 

In addition, the "whole" pressure is partitioned into terrain and 
reference parts according to 

P w -P t + Pr (A. 16) 

The hydrostatic equation is partitioned into four groups of terms. 
These are: 

Terrain, 

[ ( Ez. - 1 ) dfl w , r + dlnp r j 

v p R da da p R da 


(A. 17) 
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Reference, 

‘tg^j 

Meteorological , 

Residual, 


[ ( — -1) (4>-<t>°) + 

Pr da 


T~T° dlnp T j 
~r 


(A. 18) 


(A. 19) 


(A. 20) 


where , 

v ain Pj? 

Y " da 

Non-derivative p w and p R in (A. 17) and (A. 20) are layer mean 
pressures which must be accurately known for the partition to be 
successful. After some experimentation, it was found that, given 
the pressures at the top and the bottom of the layer, the average 
of the arithmetic mean plus twice the geometric mean, 

0.5 ( P t +P b ) +2 y[PtP~ b 
P “ 3 

yields accurate layer mean pressure. The superscript zero 
identifies observed variables. These are not subject to the 
variational operations. 

Upon specification of p R (p R = 1000, 900, 800 mb), p T is known 
through (A. 16). Therefore, (A. 17) can be solved for the terrain 
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height <p 1 . <p R is found from the level average of height after the 
removal of <f> 1 . Remaining reference variables are obtained through 
(A. 18) and the meteorological variables are found from (A. 15). The 
residual group (A. 20) exist through small modifications in <p that 
result from the variational operation. These terms are typically 
two orders of magnitude smaller than the meteorological terms. If 
these terms are represented by B, then the hydrostatic equation 
that is subject to the variational operation is, 

-i^+Y^+p-O (A. 21) 

do 


Now the pressure gradient terms of the horizontal momentum 
equations can be partitioned to separate the terrain part from the 
meteorological part that is subject to the variational operations. 
The modified nondimensional pressure gradient term is, 


d<fr r dlnp _ d<|) 
dx dx dx 


+n 


X 


(A. 22) 


where , 


dlnp T l d^ T ^ rrQ dlnp T 
dx dx w dx 


A. 3 Partition of the Local Tendencies of u and v 


Local changes in the horizontal velocity components result 
from translation of existing disturbances and development. 
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Consider that the local change in the u-component of the wind for 
a moving weather system is, 


du 

at 


_C-Vu+ — 
dt 


(A. 23) 


where c is the velocity of an advective or steering current 
(Fjortoft, 1952) usually a smoothed middle tropospheric wind. Let 
u = u 0 + u' where u 0 is the u-component of the steady part of the 
circulation and u' arises from development. Then, 


du 

dt 


-c-Vu 0 + ( 


du 7 

dt 


-C'Vu') 


The first term is the local change in u caused by translation of 
the steady part of a disturbance. The second terra is the local 
change of u from development. Note that the vertical advection of 
u is considered part of development. 

The use of the advective current throughout the troposphere is 
valid because most synoptic systems tend to maintain vertical 
structure. Any changes in vertical structure are assumed to be the 
result of development. However, the variational operations require 
that the adjustments be done on total velocity components. 
Therefore, we represent the local tendency of u by (A. 23). The 
total derivative, an approximate developmental component, is 
defined as a new dependent variable, e u = du/dt (e y = dv/dt) . 
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A. 4 The Dynamic Constraints 


Subjecting the dynamic equations (A.l) - (A. 5) to the required 
transformations yields the following constraints: For the 
horizontal momentum equations, 


R 0 [e u+ m(u-c x ) m{v-c y ) -^+* 0 d-2H] 

-a+j^av+a+j?!*) (J£+ti x )+/ u -o 


R 0 [e v +m(u-c x ) ^+m{v-c y ) ^f.+R 0 d-^] 


+ (1 +R 1 C) u+d+^JO ( +T| y ) +f v -0 


(A. 25) 


dy ° do 

W' y) 

As part of the nondimensionalization, the Coriolis parameter and 
the map scale factor have been expanded into a Taylor series. 
Thus, f = 1 + R,C and m = 1 + R,K where R 1 = 0.1. 

The continuity equation will become an integrated constraint, 


/«» ( it d °* (4 ' 4 ° > + / ff » [ 

dK 

dx 


+ R 1 K(^ + -^) -R, (ujg + vjZ) ] do - 0 


(A. 26) 


The hydrostatic and thermodynamic equations are, 


-l^+yr+p-o 

da 


(A. 27) 
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R, [« r + “<u| 1+vjZ ) *0 ( -*<1,0 
-*«!<■> .(-y r,+T) ] +*<f„--£-o 


(A. 28) 


where , 


*2 , &r 
F 1 da 


- k ^ 3 t' r ) 


is the static stability. Here F is the Froude number and Q* 
carries nondimensionalization constants. In addition, 


K- 


R 


e 



where the latter is introduced as a dependent variable. 
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Appendix B: Finite Difference Equations for the Dynamic Constraints 


The dynamic equations will be written in centered differences 
on an Arakawa D grid (Mesinger and Arakawa, 1976) . Fig. B1 shows 
the distribution of variables on the staggered grid. Anthes and 
Warner (1978) define the horizontal finite difference operators and 
the finite averaging operators as 


(B.l) 

a * = («i + l/2,j + a i-l/2,j) /2 

a (Cli,j +1 /2 + a i,j-l/ 2 ^ /2 

The i are the east-west indices, the j are the north-south indices 
as defined at the grid origin located at the lower left corner of 
the grid. In addition, the vertical differences and averages are 
defined by 


W2-°*-i/2>/Ao 

«° S ( a k*l/2 + &k-l/2) / 2 

The finite difference equations for the 
equations are, 


(B. 2) 

horizontal momentum 


M 6 —v+4> x +F 5 -0 


(B. 3) 
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M 1 ~-F 5 +R 0 {e* y +m x (u-c x ) xy u£+m x {v-c y ) u y 

+Rjo ya uZ ya ] -R 1 C*V+R l K*Q x + (1 +R 1 K X ) tl x +f u ~ 0 


M,-U+$ y +F 6 - 0 


M 2 --F 6 +R 0 [e^ +m y (u-c x ) v/+m y ( v- c y ) ^v* 
+ R 0 o x °vZ ya ] +R 1 C y u+R 1 K y 4> y + (1 +R 1 K y ) r| y +f v - 0 

The continuity equation is 

M 3 ~ / ^5 ( u x +v y ) da + o) + f Qs F i do “° 

iiiSi F n -/ [ q 2 ‘ w s +i? 1 P y ( u x + v y ) ( u ■ **/+ ] do 

The hydrostatic and thermodynamic equations are, 

M 4 “4>o + Yr+P-0 


M 5 -R 0 [e T +/n^u w T/+/n^7^r y *+d ( T^-Kq^T^) 

-Kq^to s ( T r + T* y ) ] +6o o --^--0 
F c n 


(B • 4 ) 


(B.5) 


(B.6) 


(B.7) 


(B.9) 


(B. 10) 
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The seven dynamic equations are referenced at, respectively, M 1 and 
M 6 at v, M 2 and M ? at u, M 3 at D, M 4 at T, and M g at the vertical 
velocity. 
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Appendix C: The Euler-Lagrange Equations 

The gridded fields of meteorological data to be modified are 
meshed with the dynamical equations through Sasaki's (1970a) 
variational operations. To simplify the derivations, the 
frictional terms in the horizontal momentum equations and the 
diabatic heating term in the thermodynamic equation were set to 
zero. 

Early experiments with this method found that the divergent 
part of the wind was decoupled from the adjustment with the result 
that the continuity equation was not satisfied. Attempts to 
readjust the winds through a subsidiary variational formulation 
that satisfied the continuity equation were not successful. The 
vertical velocity tended to "drift" with the result that the 
thermodynamic equation was not satisfied. 

Analysis of the problem revealed that the divergent part of 
the wind could be coupled with the variational adjustment if an 
additional constraint was satisfied. The adjusted variables must 
satisfy a particular solution of the integrated vorticity equation. 
The integrated divergence and the integrated vorticity theorem must 
vanish at the top of the model domain. This requirement is met if 
F 5 and f 6 are made dependent variables and M 3 is modified to 


AT 3 - - gr 5 (F 6x -F 5y ) A o+ (6-6 Q ) +g 5 F 7 Ao=0 


(C.l) 
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In addition. 


M 6— v+$ x +F 5 -0 


(C.2) 


Mj-U+fyy+Fs-O (C.3) 

The finite difference analog of the adjustment functional is, 

(c.4) 

i j 

The integrand, I, , is 

1 9 J 

( u - u °) 2 + n 1 ( v-v°) 2 +ti 2 (d-d °) 2 
+n 3 (<J>-<|> 0 ) 2 +ti 4 (r-r°) 2 +n 5 (<|) x -4>°) 2 
+ti 5 (<j) y -<J)p 2 +ti 6 (<j> ff -<|> 0 o ) 2 +ii 7 (e u -e2) 2 (C. 5) 

+7i 7 (e v -e °) 2 +7i 8 (e T -e £) 2 +tc 9 (F 5 -F s °) 2 

+n 9 (F 6 -F 6 °) 2 +2j>iMi 

The weights, 7T,, are Gauss' precision moduli (Whittaker and 
Robinson, 1926) . The gridded initial variables (u°, v°, a 0 , $°, T°, 
€ u °, v °, F 5 °, F 6 °) enter in a least squares formulation and receive 
according to their relative accuracies. The strong constraints to 
be satisfied exactly are introduced through the Lagrangian 
multipliers A... 

Objectively modified meteorological variables are determined 
by requiring the first variation on F to vanish. A necessary 
condition for the existence of a stationary set is that the 
functions are determined from the domain of admissible functions as 
solutions of the Euler-Lagrange equations. The variation is to be 
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carried out at every point (r,s) within the grid. Thus, upon 
setting the weights aj = bj = 1 and differentiating the integrand 
(C. 5) with respect to the arbitrary variable a r e , the Euler- 
Lagrange operator in finite differences is 


dI i,j . j 

da da 

s 1,3 




(C.6) 


Each term in I, , that contains an overbar term, that is, each term 
in Mj [ (B. 4) , (B.6), (B. 9) , (B.10), (C.l) - C.3)] produces an 
over bar term when subjected to the operations specified by (C.6). 
Multiplicate overbar terms such as ('*?) are treated having no 
overbar so that fewer grid points are required to express these 
terms in the Euler-Lagrange equations. 

The Euler-Lagrange equations resulting from the operations 
specified by (C.6) are 

T^U+Xj+F-l-O (C . 7 ) 


itiV-Xg+Fj-O 


(C.8) 


Tt 2 d+A, 3 +F 3 -0 


(C.9) 
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71 f'V 2 <J) + 7t I <|) „ „ + n 5 x (j> K 5y <j) ^ + 7t 6 „ <|) " 

- 71 3<J) + A 5jf + Xg y + A. 4 g + F 4 - 0 


7t 4 r+YA, 4 +F 8 -o 


n 7 (e u -0 +£ o X 1 -0 


^7 ® + 'R 0 ^ , 2 — ® 


Tig (e T -e?) +R 0 X 7 -o 


* 9 (^5-^5°) -A-i+Xg-Ao (g 5 X°) y -0 


ti 9 (P’g Fg ) X 2 +X 6 +Ao (g 5 A 3 ) x -0 

Variation on the Lagrange multipliers restores the 
constraints [(B.4), (B.6), (B.9), (B.10), (C.l) - C.3)]. 

The forcing functions, F 1 - F 4 contain the following 

F i°~ n i u °-AaX 3 x +R 1 C y X 2 +Rjm y X^u*+m y X 2 v x 

- [m»Tr(u=c- x ) x ] x - [ml^v-Cy) x ] y -R o (a*I?"> , 

+m y {\ 1 T*) xo }-R 1 [ (X 3 P y ) x +X 3 K^] 


(C.10) 

(C.ll) 

(C.12) 

(C . 13 ) 

(C. 14 ) 

(C. 15) 

(C . 16) 
original 


(C. 17) 
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F 2 --n 1 v 0 -/XoX 3 y -R 1 C x X 1 +Rjm x \ 1 Uy+m x X 2 y Vy 

- [mX 2 (u^) n x - [m*y\Z (v^cTy) y ] y -R 0 ) 0 ( C . 18 ) 

+M x (\ 7 T*) ya )-R 1 [ (\ 3 K xy ) y +X 3 K^] 


F 3 ~n 2 6 °+X 5 o a +R 0 X 5 ( T^ ya -K^ y T xy ) +R 2 Q (X y ° u/+Xf v/) (C. 19) 

f 4" n J*<.- 1l 5 V2 4 >0 - 7t 6 < t > ao-’t5x < l > x X -^5y < l > y y -’ l 6a4>o C>0 + *i*( *i* + * 2 y> (C.20) 

In addition, the forcing function F 8 is, 

F s — n i T°-R [ (m x u x Xj) y x + (m y v y X°) y + (SZJj f ra +Kg 3 (oXJ - ) ^ 

(C. 21) 

+Kg 4 (w S X 7 ) 

We observe that the forcing functions contain the nonlinear 
terms of their respective equations. Further, the forcing 
functions consist of terms that are either observed and therefore 
not adjusted, or are multiplied by R 0 or R 1 . These equations may 
be therefore linearized and a solution obtained through a cyclical 
method as follows. Terms multiplied by R 0 or R, are expressed with 
observed variables at the first cycle, and are expressed by 
previously adjusted variables at higher cycles. Therefore the 
forcing functions are known at each cycle. This solution method is 
valid for the latitudes and motion scales for which the Rossby 
number is less than one. 



94 


The set of equations [(B.4), (B.6), (B.9), (B.10), (C.l) - 
C . 3 ) , (C. 7) - (C. 16) ] are the linear algebraic and partial 
differential equations to be solved. Variables may be eliminated 
to reduce the number of equations to three diagnostic equations in 
vorticity, divergence, geopotential. Eliminate A 4 , A 5 , A 6 , and T 
between, respectively, (B.9), (C.10) and (C.ll); (C.8) and (C.15), 
and (C. 7) and (C.16). Next, eliminate 3 between (C.9) and (C.15) 
and (C.16). Then, A, and A 2 may be eliminated between (C.12) and 
(C. 13 ) and (C.10), (C.15) and (C.16). If M 1 and M 2 are rewritten, 
pulling out the e u and e v terms and designating the remaining terms 
as f 5 and f 6 , respectively, then e u and e v may be eliminated by 
substituting (C.12) and (C.13) into (C.15) and (C.16). Finally, 
letting D = u v + v„, the vertical velocity can be eliminated between 
(C.l) and (C.15) and (C.16). Performing the above operations 
reduces the Euler-Lagrange equation set to the following five 
equations : 

-v+4> x +F 5 -0 (C. 22) 

U+<J> y +F 6 -0 (C. 23) 


itx v + (% 9 +Jh-)F 5 - (A a) 2 [ (g|n 2 ) °D] +G, y +G 3 =0 

R I 


(C. 24) 


C-2. 
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-«1+ («9 + — I ) F 6 + (Ao) 2 [ (g5 2 ^ 2 ) X- G 1X- G 2"° 

Ro 


(C. 25) 


7tf'V 2 <|> + 7l?<l> <y< , + 7l 5jc <j)^+Jl 5y <|)^+lI 6o 4)" 
-71 3 4>+ (t^v) x - (n 1 u) y +G 4 -0 


(C. 26) 


where the forcing functions, G 1 - G 4 are given by: 


Gi-Ao [q 5 F°+q 5 n 2 6°-q 5 nlF7bo] 


G 2 -* 9 F 6 °+^(f'+R ( fil)+F x 

Rl 


G 2 —n 9 F 5 °-^-(f 5 +R 0 e° u )+F 2 


Gt-F^-F, „+ ( Z4- + 2±E1) o+ f 4 


'4 x 2x * iy • y2 y 


We are now in a position to substitute (C.22) and (C.23) into 
(C. 24) and (C. 25) to eliminate F 5 and F 6 . We make note that the 
substitution generates the following combination of precision 
modulus weights, 

Further, we note that all of these precision moduli vary 
horizontally with horizontal variations in Tt y . Thus, if, 
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Ro Ro 

(X/Y/ a ) =7r i ( a ) f (X/Y) f and the horizontal variations of n 7 and n 9 
also vary as f(x,y), then by dividing all precision moduli by 
f( x /Y)/ the horizontal variations of 7r 10 and 7r n may be removed 
without changing the relative relationships between the weights. 
With these modifications, the Euler-Lagrange equations (C.24) and 
(C. 25) may be combined to form a divergence equation, 

V 2 [ (Ao) 2 g|jt 2 I?] - it 10 Z>= V 2 + G 2x + G 3y (C.27) 

The vorticity formed from (C.24) and (C.25) is, 

(C.28) 


Substitution of the vorticity between (C.26) and (C.28) leaves a 
diagnostic equation in geopotential, 


(^f'+ni^il ) ^<1)+ (Ttg + ^f ) <j> 0O + it 5jc < |)^+ii 5y <|)^ (ti 6o +-^|- ) <j>"-7C 3 <|) 


10 


n 


-G A --±(G 2y -G 3x ) 


(C.29) 


n 


10 


Equations (C.27) - (C.29) form the three diagnostic equations that 
must be solved for a successful variational adjustment. All terms 
to the right of the equal sign are forcing functions that contain 
either unadjusted initial variables and/or variables that have been 
adjusted at the last iteration. (C.29) is solved first to get the 
geopotential height. Then the divergence and vorticity are 



obtained through (C.27) and (C.28). 
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Appendix D: Boundary Conditions 


The variational theory specifies natural boundary conditions 
that are consistent with the Euler-Lagrange equations. If it is 
assumed that there are no adjustments in the data along the 
boundaries, then the boundary conditions may be specified. In the 
latter case, the Lagrange multipliers. A.,., are zero at the 
boundaries and the initial unadjusted values are used for the 
boundary conditions. 

Initially, the Euler-Lagrange equations were solved with 
specified boundary conditions. These boundary conditions forced 
high frequency waves into the solutions for the velocity components 
near the boundaries. Divergences calculated from these velocity 
components gave large erroneous vertical velocities. We therefore 
returned to the natural boundary conditions. 

The Euler-Lagrange operator for natural boundary conditions 
is. 



(D.l) 


Performing the operation specified by (D.l) yields the following 
expressions for the boundary conditions on $ 


n 5 (<t> x -<|>x) +A 5 +i? 1 i?A. 1 =0 


(D.2) 
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^5 ^y) Xg +^3^X3“ 0 


(D.3) 


* 6 (*.-*.°>+V0 (D * 4) 

The terms multiplied by R 1 come from the constraints, M 1 and M 2 . 
These equations can be solved for the § boundary conditions subject 
to substitutions for the X f through the Euler-Lagrange equations 
(22) -(35) in the text. The lateral boundary conditions for the x- 
and y-boundaries are, respectively, 


it 


(*10*5 + *13 ^ *io*5^x + ~ — •^2 + *9*12'^5 + ^12^ ® ) 


Tt, 


+ ^I (i+n 12 ) (-^-+eS) -o 


(D. 5) 


It 


(it 10 7t 5 +it 13 )<l> y -7t 10 ii 5 <j)°-— ^■F 1 +it 9 7i 12 F 6 0 -n 12 Ao((?X3) x 


71 


+ 4 1 (1+* 12 > (-r +e ^ “° 


(D.6) 


where , 


*io- 


it? 

*9 + *l + ^ 

It. 


71? 


R.K 7t 7 

71., “1 ~ — , 

12 *1 




k 13 11 12 


10 1 
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Several observations may be made with regard to (D.5) and 
(D . 6) . 

(1) F 1 , F 2 , 1 3 , f 5 , and f 6 all contain terms that are updated at 
each cycle. Thus it is possible to update the boundary 
conditions as the interior fields are being adjusted. 

(2) These forcing functions contain nonlinear terms that 
cannot be calculated at the boundaries unless derivatives are 
extrapolated across the boundaries. Therefore, the boundary 
equations may be simplified by setting A >1 = A 2 = X 3 = 0 at the 
boundaries. It follows therefore, that 

F x — itjU 0 , F 2 --n 1 v°, F 2 ~-n 2 6° 

(3) From (22) and (29) , 

Jl, 

QAo (D+F 7 ) +— 1+6 °-d 0 -0 (D . 7 ) 

n 2 

Given that it is the gradient of A 3 that appears in (D.5) and 
(D. 6) it follows from (D.7) that gradient of the divergence 
must be specified, or in other words, the divergence must be 
specified along at least two boundary grid rows or columns in 
order that the gradient of 1 3 vanish in the $ boundary 
equations. 

(4) 7r 7 is at least two orders of magnitude smaller than the 
remaining precision moduli. Neglecting n 7 leads to the 
following simplifications, 

The equations for the lateral boundary conditions on $ are thus, 
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Ttg+TCi 

*10- — / 


'll n 13 " 9 ' 12 


*1*9 


7C ^ TC c 


( TC + ) 4> -7I 5 4> ° 3^-L- ( U°-F 5 ° ) =0 

5 7l 1 + 7l 9 - ■ - 


7li + 7l 9 


(D.8) 


(it 5 + 7tlTC — )<{> -Tt 5 $°+ 

5 7^ + Ttg y 5 y 7^+ltg 


o. *1^ ( u o +j p 6 °) , 0 


The boundary conditions for u and v may be found by solving 
the same set of equations used for finding the $ boundary 
conditions but for u and v. The results are, 


(715 + 11! 
(715 + 71! 


7t 5 + 7I 9 
7t a 


7t 5 + 7t 9 


) V-7T 5 (4>°-F 5 °) — 7t ! IhUh v°-0 

1t 9 

) u+7t 5 ($y+F<°) -Tt! rc 5+ * 9 U °-0 

Ttn 


(D.9) 
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Table 1 

Nondimensional standard errors of observation for wind, height, 
and temperature and RMS errors for other adjustable meteorological 
variables. 


VARIABLE 

Model Pressure Mean 

Level (mb) u 20 u 40 $ A$/Ax A$/Aa Temp a e u 










0.00 


10 

100 

0.45 

0.23 

0.25 

0.71 











3.68 

0.59 

2.13 

6.98 

9 

200 

0.45 

0.23 

0.20 

0.56 











3.21 

0.88 

1.88 

6.98 

8 

300 

0.42 

0.21 

0.18 

0.51 











2.28 

0.88 

1.64 

6.51 

7 

400 

0.36 

0.18 

0.15 

0.42 











1.53 

0.76 

1.43 

5.58 

6 

500 

0.32 

0.16 

0.12 

0.33 











0.97 

0.59 

1.24 

4.65 

5 

600 

0.30 

0.15 

0.09 

0.26 











0.61 

0.44 

1.04 

4.34 

4 

700 

0.28 

0.14 

0.08 

0.22 











0.53 

0.44 

0.84 

3.72 

3 

800 

0.24 

0.12 

0.07 

0.20 











0.47 

0.44 

0.64 

3.26 

2 

900 

0.21 

0.11 

0.06 

0.18 











0.42 

0.44 

0.44 

3.10 

1 

1000 

0.20 

0.10 

0.06 

0.17 
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Table 2 


Nondimensional precision modulus weights for variational objective 
analysis. 








VARIABLE 




Model 

Pressure 




Mean 




Level 

(mb) 

u 20 

$ 

At/Ax At/A a 

Temp 

a 

e u 









100 



10 

100 

2.5 

8.0 

1.0 

0.04 

1.4 

10 

0.01 


9 

200 

2.5 

12.5 

1.6 

0.05 

0.6 

0.14 

0.01 


8 

300 

2.8 

15.4 

1.9 

0.10 

0.6 

0.19 

0.01 


7 

400 

3.9 

22.2 

2.8 

0.21 

0.9 

0.24 

0.02 


6 

500 

4.9 

34.7 

4.6 

0.53 

1.4 

0.33 

0.02 


5 

600 

5.6 

61.7 

7.4 

1.34 

2.6 

0.46 

0.03 

• 

4 

700 

6.4 

78.1 

10.3 

1.78 

2.6 

0.71 

0.04 


3 

800 

8.7 

102.0 

12.5 

2.26 

2.6 

1.22 

0.05 


2 

900 

11.3 

138.9 

15.4 

2.83 

2.6 

2.58 

0.05 


1 

1000 

12.5 

138.9 

17.3 
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Table 3. Correlation coefficients for a 216-point 
subset of initial (i) and variational (v) u, v, 
and T 3-h forward tendencies at 0000 UTC compared 


with observed 3-h tendencies centered at 0130 UTC. 


p 

lev 

u. 

U v 

_ Vj 

V 

T i 

T y 

200 

-0.34 

V 

-0.08 

-0.27 

V 

-0.25 

0.17 

0.07 

300 

-0.10 

-0.24 

0.43 

0.10 

-0.36 

0.17 

400 

-0.24 

0.12 

0.53 

0.35 

0.24 

0.59 

500 

-0.26 

0.31 

0.43 

0.71 

0.75 

0.65 

600 

0.36 

0.56 

0.01 

0.35 

0.42 

0.75 

700 

0.66 

0.71 

0.15 

0.61 

0.55 

0.72 

800 

0.55 

0.59 

0.54 

0.79 

0.48 

0.17 

900 

0.65 

0.60 

0.31 

0.37 

0.25 

0.22 


Table 4. Same as Table l but with 0000 UTC 
3-h forward tendencies shifted by weather system 
translation to approximate 0130 UTC observed 
tendencies. 

P 

lev 

u i 

U v 

V i 

V u 

T i 


200 

-0.35 

V 

-0.06 

-0.31 • 

v 

-0.25 

1 

0.12 

0.02 

300 

0.01 

-0.04 

0.56 

0.23 

-0.31 

0.14 

400 

-0.27 

0.20 

0.45 

0.30 

0.35 

0.67 

500 

-0.03 

0.57 

0.54 

0.73 

0.83 

0.64 

600 

0.56 

0.69 

0.02 

0.45 

0.57 

0.76 

700 

0.79 

0.83 

0.22 

0.73 

0.66 

0.71 

800 

0.72 

0.75 

0.60 

0.87 

0.55 

0.23 

lll'IMW 

0.82 

0.78 

0.35 

0.48 

0.32 

0.49 
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FIGURE CAPTIONS 

Fig. l. The distribution of rawinsonde stations over the analysis 
grid (solid rectangle) , evaluation grid (large dashed 
rectangle) , and SESAME I network (small dashed rectangle) . 

Fig. 2. Residual reduction as a function of cycle for the u- 

component (left panel) and v-component (right panel) dynamic 
constraints. 

Fig. 3. Residual reduction as a function of cycle for the 
integrated continuity equation (left panel) , the hydrostatic 
equation (middle panel) , and the thermodynamic equation (right 
panel) . 

Fig. 4. RMS differences between unadjusted (adjusted) fields and 
observations after removal of standard observation error 
(solid lines) and means of differences between unadjusted 
(adjusted) fields and observations (dashed lines) for a) 
heights, b) temperatures, c) u-comp, and d) v-comp. 

Fig. 5. Heights and wind vectors at 800 mb, 500 mb, and 300 mb for 
a) unadjusted and b) adjusted fields. 

Fig. 6. Differences between adjusted and unadjusted heights and 
vector winds at 800 mb, 500 mb, and 300 mb. 

Fig. 7. Same as Fig. 5 but for temperature. 

Fig. 8. Relative vorticities at 500 mb, a) unadjusted and b) 
adjusted. 

Fig. 9. a) unadjusted, b) adjusted vertical velocities (cm sec' 1 ) 
at 500 mb. Precipitation areas are stippled. 

Fig. 10. u-component tendencies for 800 mb (left panels) and 500 mb 
(right panels) for a) observed, b) unadjusted, and c) adjusted 
fields in m sec* 1 3-hr' 1 . 

Fig. ll. Same as Fig. 10 but for the v-component. 

Fig. 12. u-component tendencies (left panels) and v-component 

tendencies (right panels) at 300 mb for a) observed, b) 
unadjusted, and c) adjusted fields in m sec' 1 3-hr' 1 . 




Fig. 1. The distribution of rawinsonde stations over the analysis 
grid (solid rectangle), evaluation grid (large dashed 
rectangle), and SESAME I network (small dashed rectangle). 


RESIDUAL 




Fig. 2. Residual reduction as a function of cycle for the u- 

component (left panel) and v-component (right panel) dynamic 
constraints . 





Fig. 3. Residual reduction as a function of cycle for the integrated 
continuity equation (left panel), the hydrostatic equation 
(middle panel), and the thermodynamic equation (right panel). 
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Fig. A. RMS differences between unadjusted (adjusted) fields and 

observations after removal of standard observation error (solid 
lines) and means of differences between unad justed (adjusted) 
fields and observations (dashed lines) for a) heights, b) 
temperatures, c) u-comp, and d) v-comp. 
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Fig. 5. Heights and wind vectors at 800 mb, 500 mb, and 300 mb for 
a) unadjusted and b) adjusted fields. 
















a 


b 


300 mb 




500 mb 



Fig 
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Same as Fig. 5 but for temperature 










Fig. 8. Relative vorticities at 500 mb, a) unadjusted and b) 
adjusted. 







800 mb 


500 mb 



Fig. 10. u-component tendencies for 800 mb (left panels) and 500 mb 
(right panels) for a) observed, b) unadjusted, and c) adjusted 
fields in m see"* 3-hr - *. 
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Fig. 11. Same as Fig. 10 but for the v-component. 










Fig. 12. u-component tendencies (left panels) and v-component 
tendencies (right panels) at 300 mb for a) observed, b) 
unadjusted, and c) adjusted fields in m sec 3-hr 
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1. Introduction 

The MODEL III variational data assimilation model is the third 

of four general assimilation models designed to blend weather data 

measured from space-based platforms into the meteorological data 

mainstream in a way that maximizes the information content of the 

satellite data. Because there are many different observation 

locations and there are many instruments with different measurement 

error characteristics, it is also necessary to require that the 

blending be done to maximize the information content of the data 

and simultaneously to retain a dynamically consistent and 

reasonably accurate description of the state of the atmosphere. 

This is ideally a variational problem for which the data receive 

relative weights that are inversely proportional to measurement 

error and are adjusted to satisfy a set of dynamical equations that 

govern atmospheric processes.?! Because of the complexity of this 

i . . 

type of variational problem, we have divided the problem into four 
variational models of increasing complexity. The first, MODEL I, 
includes as dynamical constraints the two horizontal momentum 
equations, the hydrostatic equation, and an integrated continuity 
equation. The second, MODEL II includes as dynamical constraints, 
the equations of MODEL I plus the thermodynamic equation for a dry 
atmosphere. MODEL III includes the equations of MODEL II plus the 
radiative transfer equation. 
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The advantage of MODEL III over the previous two models is 
that radiance, the atmospheric variable measured by satellite, 
becomes a dependent variable. In the previous versions, mean layer 
temperatures that had been retrieved from the radiances by some 
method, were included in the assimilation by substituting them in 
place of the rawinsonde temperatures. Now both rawinsonde 
temperatures and satellite radiances are included independently in 
the assimilation. 

Our approach to the development of MODEL III has been to 
divide the problem into three steps of increasing complexity. 
Chapter IV deals with the first step, a variational version of the 
classical temperature retrieval problem that includes just the 
radiative transfer equation as a constraint. The radiances for each 
of the four TOVS MSU microwave channels are dependent variables. 
These plus temperature constitute a set a five adjustable 
variables. Each radiance is related to the temperature through its 
radiative transfer equation. There are therefore four dynamic 
constraints in this first variational problem. 

Chapter V summarizes the second step which combines the four 
radiative transfer equations of the first step with the equations 
for a geostrophic and hydrostatic atmosphere. This step is intended 
to bring radiance into a three-dimensional balance with wind, 
height, and temperature. The use of the geostrophic approximation 
in place of the full set of primitive equations allows for an 
easier evaluation of how the inclusion of the radiative transfer 
equation increases the complexity of the variational equations. 
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The third and final step includes the four radiative transfer 
equations with the fully nonlinear set of primitive equations, ie., 
MODEL III. 

2 . A Variational Retrieval Algorithm 

The radiative transfer equation is the only variational 
constraint. It takes the form 

OQ 

B-B 0 w 0 -fw'Tdz - 0 (1) 

o 

where B is the brightness temperature as computed from radiance 
measured at the satellite and T is the mean layer temperature of an 
incremental depth of the atmosphere, dz. The weight, w 0 , is the 
transmittance of the total atmosphere from the surface (where the 
surface brightness temperature, B 0 , is measured) to the space-based 
observation platform. The weights, w', are proportional to the 
transmittance from some level within the atmosphere to the 
satellite. In order to make the variational derivations from (1) 
compatible with the larger set of variational equations in MODEL I 
and MODEL II, we will make the following modifications in (l). 
First, the brightness temperature is replaced by the skin 
temperature, T 0 , and the weight, w 0 , will become a skin level 
surface weight. Second, (1) is converted from the z to the sigma 
vertical coordinate. In this conversion, 
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J w'Tdz- J t ^T[f(T) ] da- J wTda 


( 2 ) 


Now f (T) , a small conversion term that results from the changeover 
to sigma coordinates, will be combined with the weights and not 
subjected to variation. This approach avoids complicated nonlinear 
equations that will otherwise arise through the variational 
formations. The f(T) and the weak temperature dependence in the 
weights will not be held constant however. At each step of a 
converging iterative process, the small temperature dependencies 
will be updated with adjusted temperatures. With these 
modifications, (1) becomes, 

oo 

B-fwTdo- 0 (3) 

o 

The next step is to bring (3) into dimensional compatibility with 
the more general variational models. Let, 

T-Qt'-SHt' (4) 

R 


and 


T'-T r + 


_F_ 


rj*H 


(5) 
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so that. 


T~-^ (T r + — T") 


R 


R 


( 6 ) 


Here g is gravity, H=10 km is a reference height, R is the 
universal gas constant, F is the Froude number, and R 0 is the 
Rossby number. The subscript R refers to a reference atmosphere 
and the notation " refers to departures from the reference 
atmosphere. Substitution of (6) in place of T in (3) gives, 

00 oo 

B--&H [ fwT R do+— fwr"do] =0 (7) 

R J o R °o 


Further, we partition B = B R + B m and define 



( 8 ) 


It follows then, that 


gH_F. 

R R 


( 9 ) 


Finally, upon suppression of the double primes, the radiative 
transfer equation becomes, 
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K 


B- E 

Jc-1 


( 10 ) 


Now there are four TOVS microwave channels each with an 
independent measurement of the brightness temperature. Let Bj be 
the brightness temperature perturbation for the jth channel. The J 
constraining equations are, 

K 

m r B rE w * v T *’° (11) 

Jc-l 


The functional to be minimized is 

f Ida 


( 12 ) 


where 


K 


j 


f-E * * ( T k - T k°) it , 

k-1 j - 1 


*E if 


t m. 


J-l 


(13) 


Performing the variations upon T and B as shown by Achtemeier, et 
al. (1986) yields the following Euler-Lagrange equations, 

6T: n k (T k -T k °)-J2w kj X r 0 


(14) 
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for each k and, 


6 Bi Kj (Bj-Bf) +Xj-0 


( 15 ) 


for each j . Variation upon the J Lagrange multipliers restore the 
original constraints (11) . These equations are linear and may be 
easily reduced to one diagnostic equation in temperature. First, 
eliminate reference to the Lagrangian multipliers by substituting 
(15) into (14) . Then substituting for Bj gives the adjusted 
temperature as a function of weight functions and observed 
variables, 

\ 

n * T * + £ £ n 3 W U W ii T r F k- 0 <16) 


for each k. Here 

**- w * r * + £ * jBf 

j-i 


( 17 ) 


Equation (16) can be easily solved with a standard matrix 
inversion package to retrieve the variationally adjusted 
temperature profile. At most two cycles with the weight functions 
updated with adjusted temperatures are required for convergence to 
a final adjusted temperature. 
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4 . Results 

In order to properly interpret the results of the example of 
variational adjustment with the radiative transfer equation, one 
must be aware that three sets of weights appear in (16) . The 
weights, w^. , are the transmittance weights for the ith level and 
the jth microwave channel. They are not subject to the variational 
adjustment and remain unchanged with the exception of minor 
adjustments for temperature sensitivity. The variational weights, 
TTj and 7T k , carry the relative importance of the jth microwave 
channel and the temperature at the kth level. It is the choice of 
the variational weights that are important in interpreting the 
results. 

Consider a temperature profile that is to be retrieved from 
MSU brightness temperatures. It is to be made halfway between two 
rawinsonde sites. The rawinsonde soundings are given by A and B in 
Fig. 1. Sounding A is cold up to the tropopause (about 220 mb) and 
then it becomes isothermal up to 60 mb. Sounding B is warm from 
the surface to 170 mb and then becomes colder than A in the layer 
from 170 mb to 60 mb. Its tropopause is located at 100 mb. 

The first guess or "observed" sounding that will enter into 
the temperature part of the variational analysis is the mean of A 
and B. It is given by M in Fig. 1. Now suppose that the true 
sounding is given by T. Note that M=T from the surface to 230 mb 
and from 50 mb to the top. 
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Next, the brightness temperatures, Bj, were calculated from 
(10) using the true temperature sounding. Thus the Bj° that enter 
(17) are true and the T k ° are approximate. However, only the 
temperatures between 100-230 mb need adjustment. The observational 
error for the temperature was 0.7 K and the weight accorded to the 
temperature was, 

n k ±-- 1.0 (18) 


Fig. 2 shows the results of three retrievals between 500 mb 
and the top. The dashed line is the difference M-T between the 
true and first guess temperature soundings. The other curves are 
the differences between the adjusted and the true temperatures for 
7Tj that ranged in values from 10 to 100 to 1000. Note that the 
weights for the four MSU channels and hence the brightness 
temperatures were always equal. 

Fig. 2 shows that increasing the brightness temperature 
weights progressively reduced the differences between the adjusted 
and true temperature soundings but by only 2.5 K. However, the 
retrievals also spread the adjustments throughout the depth of the 
sounding. Therefore, improvements where the M-T residuals were 
nonzero were offset by degraded temperatures throughout the 
reminder of the sounding - the errors being almost 2 K at 250 mb 
with lesser error elsewhere. 

A more extensive analysis of the behavior of (16) found that 
the retrievals were sensitive to the vertical distribution of the 
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weights for the temperature hence the errors of observation for the 
temperature. If there existed some independent observations that 
could be used to estimate the accuracy of the first guess 
temperature as a function of height, then the retrievals could be 
focused into those locations where the M-T residuals were greatest. 
Consider possible accuracy functions given in Fig. 3. The 
effective temperature error at 150 mb is doubled by f(l) and is 
tripled by f ( 2 ) . Therefore, the weights accorded to the 
temperature there are decreased by a factor of four for f ( 1 ) and a 
factor of nine for f ( 2 ) . 

Fig. 4 shows the residuals between the adjusted and true 
temperature profiles for the three retrievals when the accuracy 
function f(l) was applied to the temperature weights. The initial 
residual has been reduced by approximately 6 K. Fig. 5 shows the 
results for f ( 2 ) . Additional reductions in the residuals over f(l) 
results were found between 150 and 100 mb. Fig. 6 summarizes the 
resulting temperature soundings for f(0), f(l), and f(2) if the 
weights for the brightness temperatures were n. = 1000. The 
improvement of f ( 2 ) over f(l) is apparent between 150 and 100 mb 
but elsewhere the differences between the two retrievals are only 
a few tenths of a degree. This suggests that it is the shape of 
the accuracy function, not the magnitude, that determines where the 
variational adjustment will be focused. 

Fig. 7 shows part of the temperature soundings T and M between 
250 mb and 50 mb. The curve identified by VI is the sounding that 
was obtained with the conditions that the weights for the first 
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guess temperature were constant with height. The sounding V2 
results from the application of f ( 2 ) to the temperature weights. 

The first step in the variational analysis of the radiative 
transfer equation succeeded in producing a variational algorithm 
that could be used to retrieve temperature from the four MSU 
channel brightness temperatures given a first guess temperature 
sounding. The results showed that the variational retrievals were 
subject to the same limitations as are retrievals by other methods, 
inability to accurately resolve temperatures near the tropopause 
spreads error though the whole retrieved sounding, unless some 
temperature accuracy function is employed to focus the retrieval. 
The identification of a data set that could be used for a 
temperature accuracy function and the derivation of the same is 
beyond the scope of this study. 
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FIGURE CAPTIONS 


Figure l. Two typical temperature soundings A and B; the mean of 
A and B, sounding M; and true temperature sounding T used for 
sensitivity studies of variational temperature retrievals. 

Figure 2. Dashed line: differences between the mean or first guess 
temperature sounding and true sounding. Solid lines: 
differences between variational temperature retrievals and 
true temperature sounding for the following choices of 
brightness temperature weights; sounding 1 (10) , sounding 2 
(100), sounding 3 (1000). 


Figure 3 . 
Figure 4. 
Figure 5. 


Curves for hypothesized temperature accuracy functions. 
Same as Fig. 2 but for f(l). 

Same as Fig. 4 but for f(2). 


Figure 6. Differences between first guess and true temperature 

(dashed line) and variational temperature retrievals and true 
temperature for brightness temperature weights equal to 1000 
for f (0) , f (1) , and f (2) . 


Figure 7. Parts of temperature soundings T and M between 250 mb 
and 50 mb. Sounding VI is temperature retrieval with f (0) and 
sounding V2 is temperature retrieval with f(2). 




Figure l. Two typical temperature soundings A and B; the mean of 
A and B, sounding M; and true temperature sounding T used for 
sensitivity studies of variational temperature retrievals. 
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Figure 2. Dashed line: differences between the mean or first guess 
temperature sounding and true sounding. Solid lines: 
differences between variational temperature retrievals and 
true temperature sounding for the following choices of 
brightness temperature weights; sounding 1 (10), sounding 2 
(100), sounding 3 (1000). 
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Figure 5. Same as Fig. 4 but for f(2). 




Figure 6. Differences between first guess and true temperature 

(dashed line) and variational temperature retrievals and true 
temperature for brightness temperature weights equal to 1000 
for f (0) , f ( 1) , and f (2) . 
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Figure 7 . Parts of temperature soundings T and M between 250 mb 
and 50 mb. Sounding VI is temperature retrieval with f (0) and 
sounding V2 is temperature retrieval with f(2). 
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1 . Introduction 

The approach to the development of MODEL III has been to 
divide the problem into three steps of increasing complexity. In 
Chapter IV we successfully developed a variational algorithm for 
the classical temperature retrieval problem that includes just the 
radiative transfer equation as a constraint. The radiances for each 
of the four TOVS MSU microwave channels were dependent variables. 

Chapter V summarizes the second step which combines the four 
radiative transfer equations of the first step with the equations 
for a geostrophic and hydrostatic atmosphere. This step is intended 
to bring radiance into a three-dimensional balance with wind, 
height, and temperature. The use of the geostrophic approximation 
in place of the full set of primitive equations allows for an 
easier evaluation of how the inclusion of the radiative transfer 
equation increases the complexity of the variational equations. 

It should be noted that the variational method is a powerful 
mathematical tool and a powerful method for diagnosing the physical 
role of the observations in the adjustment. We developed seven 
different variational formulations for the geostrophic, hydrostatic 
and radiative transfer equations. The first derivation was too 
complex to yield solutions that were physically meaningful. For 
the remaining six derivations, the variational method gave the same 
physical interpretation - the observed brightness temperatures 
could provide no meaningful input into a geostrophic, hydrostatic 
balance - at least through the problem-solving methodology employed 
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in these studies. It would be axiomatic therefore, that the 
brightness temperatures could provide no meaningful input into a 
variational assimilation with the primitive equations. 

During the writing of this chapter, the equations were 
reviewed and a conceptual error regarding one of the Lagrange 
multipliers was discovered. 

In the following section, the variational methodology is 
presented and the Euler-Lagrange equations rederived for the 
geostrophic, hydrostatic and radiative transfer equations. Then 
the equations are reduced in number through elimination of 
variables to produce a single equation for the geopotential height. 
It is shown that the single equation is too difficult to solve but 
that a three equation set can be solved iteratively. It is also 
shown that space-based thermodynamic data can be assimilated into 
the meteorological data mainstream and that none of the 
difficulties associated with traditional temperature retrievals 
will be encountered. 

2. A Variational Assimilation Theory for the Geostrophic, 

Hydrostatic and Radiative Transfer Equations 

The variational formalism will be derived for the four 
radiative transfer equations in integral form. Let the dynamical 
constraints be, 


^-(JW^+p-o 


( 2 ) 
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J J o° 

(Bj-fwjTdo-0) ( 1 ) 

>i j-i o 

m 3 -v-<b x -F s -0 (3) 

/n 4 =u+<J> y +F 6 =0 (4) 

For additional simplification, set the terrain correction term 6=0. 
The forcing functions F 5 and F 6 (see Chapter II) are simplified 
through setting R 0 = 0. 

The integrand of the functional to be minimized is, 

I-k-l ( u-u °) 2 +ti 1 (v-v°) 2 +tc 2 (T-T°) 2 +ti 3 (4>-<t>°) 2 

j j ( 5 ) 

+ £ ( Bj-Bf ) 2 +2j^ A, xjMxj +2k 2 m 2 +2X 3 m 3 +2k 4 m i 

Pi Pi 


where the n i are the relative weights accorded to the observations. 

Performing the variations for the eight dependent variables, 
u, v, §, T, and Bj (j=l,4), yields the following Euler-Lagrange 
equations, 


6 u: n x (u-u °) +X 4 -0 


( 6 ) 
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dv: n 1 (v-v°) + A 3 -0 (7) 

6<(>: n 3 (4>-<J>°) -X 2a +X 3x -X 4y -0 (8) 

6r: % 2 (T-T o )+Yk 2 -^X 1J J^w j da-0 (9) 

85: g [n 4J (flj-fl/)+X w -0] (10) 


These eight equations plus the seven original constraints 
constitute a set of 15 algebraic and linear partial differential 
equations to be solved. The number of equations may be reduced 
through the elimination of variables. There results a single 
diagnostic equation with geopotential height as the dependent 
variable. We develop a diagnostic equation for the geostrophic, 
hydrostatic adjustment first and then include the contribution from 
the radiative transfer equations. Two Lagrange multipliers are 
eliminated by combining (6), (7), and (8). Then, forming the 
vorticity from (3) and (4) and combining with (8) gives. 


7i 1 V 2 <J)-7i 3 <|)+A 20 -7i 1 (v°-u°) +7t 3 4>°+n 1 (F 5x +F 6y ) =0 


( 11 ) 
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Reducing the thermodynamic variables is done as follows. 
Divide (9) by y and operate by a. Eliminate brightness temperature 
between (1) and (10). There results two equations, 

l[^(r-ni4 !t -£Vj-o (12) 


where , 


e 


a 

da 



( 13 ) 


and, 



Tda-Bf) +X 1J7 - 


-0] 


( 14 ) 


Combining (12) with (14) and substituting (2) gives, 


T y 2 f=i { V 

n. 


( 15 ) 


» 3 - 1 


Eliminating the Lagrangian multiplier between (11) and (15) yields 
a diagnostic equation in the geopotential height. 
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*1^ + “T ‘ ♦o<» + (-^7 £ >o-*3 < I >+ ]£ e J K 4jf -^&-do+F= 0 (16) 

Y 2 Y J=i J 0 Y 


where , 


F-n 1 (v°-u y °) +* 3 <t>°+ (-^r°) 0 +g e (F 5x +F 6y ) 


(17) 


Much effort was spent programming for (16) . The resulting 
solution was not considered satisfactory. Given the complex 
coefficient structures and the delicate convergence criteria, much 
additional effort was expended through six subsequent derivations 
to express the variational formalism in forms easier to understand 
and easier to solve. These efforts eventually led away from a 
direct inclusion of the radiative transfer equation in a 
geostrophic, hydrostatic atmosphere. 

In retrospect, it seems that the solution could have been more 
easily obtained if the 15 equation set was reduced to the following 
3 equation set: 

7i 1 V 2 <t>-n 3 <| > --X 20 +n 1 (v°-u y ) (F 5jf +F 6y ) (18) 


T—- ^(4>„+P) 


( 19 ) 
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j. oa oo 

\ 2 —— (T-T°) -J2 f Wjda ( f WjTda-Bf) ] (20) 

Y j-i Y o o 

These equations can be solved iteratively by first setting X 2 =0 
and solving (18) for $. Then (19) is solved for T and the 
temperature substituted into (20) to derive an updated X 2 . Then 
the updated values are entered into (18) and the cycle repeated 
until a satisfactory level of convergence is attained. 


3 . Results 

There are three important points to consider regarding (18) - 

( 20 ) . 

a) NO RETRIEVAL OF TEMPERATURE IS REQUIRED TO BLEND SATELLITE 
OBSERVED BRIGHTNESS TEMPERATURES INTO THE METEOROLOGICAL DATA 
MAINSTREAM. This means that none of the problems associated 
with temperature retrievals will be encountered. The integral 
term appears on the right hand side of (20) not as a term to 
be solved. This is analogous to solving the radiative 
transfer equation for the brightness temperature - a very easy 
exercise. This single finding may make it worth while to 
pursue the formal variational approach to assimilation of 
microwave channel data especially if higher resolution 
radiance data becomes available in the future. 
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b) There must be observations of geopotential height or winds or 
both of equivalent accuracy with the satellite measurements in 
order for the (18) - (20) to work. Accurate observations of 
temperature apart from space-based measurements are not 
necessary. The caveat is that geopotential height must be 
known at the boundaries of the domain in order to obtain a 
solution for (18) . Lateral boundaries would vanish for the 
equations written on the sphere and the top boundary 
conditions can be removed to the top of the atmosphere or to 
some level where model predictions or climatology give 
satisfactory estimates. 

c) It is highly probable that (18) - (20) converge to a solution. 
The same equation set with the absence of the second term of 
(20) (the radiative transfer equation contribution) is known 
to converge. The second term of (20) is an integral term 
which should further stabilize the solution. 

The main goal of the variational assimilation project was to 
blend satellite-derived thermodynamic data into the meteorological 
data mainstream in a dynamically consistent way. The classical 
variational calculus method used to achieve that goal typically 
yields sets of complicated equations that require innovative 
methods for solution and also involve immense programming efforts. 
Therefore the effort was broken down into several simpler models 
that could be solved. 
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The attempt to reduce the equation set from 15 equations to 
one diagnostic equation in geopotential height resulted in equation 
(16) . After an extensive programming effort, a satisfactory 
solution was not obtained. I was unable to devise a scheme that 
could determine whether the problems were mathematical or 
pr ©grammatical. During the six other efforts to derive a more 

tractable diagnostic equation a conceptual error was made, namely, 
was treated as a variable that could be differentiated with 
respect to a. The observed brightness temperature dropped out of 
the equations. This led to the conclusion that the satellite data 
could not be successfully included in a classical variational 
assimilation. 

With the discovery of this error during the writing of Chapter 
5 of this final report, that conclusion is no longer valid. It 
appears, instead, that the satellite data can be successfully 
incorporated into a variational assimilation and that the blending 
can be done without any of the problems typically encountered with 
temperature retrievals. 
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ABSTRACT 

There has been a long-standing concept by those who use 
successive corrections objective analysis that the way to obtain 
the most accurate objective analysis is first, to analyze for the 
long wavelengths and then to build in the details of the shorter 
wavelengths by successively decreasing the influence of the more 
distant observations upon the interpolated values. Using the 
Barnes method, we compared the filter characteristics for families 
of response curves that pass through a common point at a reference 
wavelength. It was found that the filter cutoff is a maximum if 
the filter parameters that determine the influence of observations 
are unchanged for both the initial and corrections passes. This 
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information was used to define and test the following hypothesis. 
If accuracy is defined by how well the method retains desired 
wavelengths and removes undesired wavelengths, then the Barnes 
method gives the most accurate analyses if the filter parameters on 
the initial and corrections passes are the same. This hypothesis 
does not follow the usual conceptual approach to successive 
corrections analysis. 

Theoretical filter response characteristics of the Barnes 
method were compared for filter parameters set to retrieve the long 
wavelengths and then build in the short wavelengths with the method 
for filter parameters set to retrieve the short wavelengths and 
then build in the long wavelengths. The theoretical results and 
results from analyses of regularly spaced data show that the 
customary method of first analyzing for the long wavelengths and 
building in short wavelengths is not necessary for the single 
correction pass version of the Barnes method. Use of the same 
filter parameters for initial and corrections passes improved the 
analyses from a fraction of a percent for long wavelengths to about 
ten percent for short but resolvable wavelengths. 

However, the more sparsely and irregularly distributed the 
data, the less the results are in accord with the predictions of 
theory. Use of the same filter parameters gave better overall fit 
to the wavelengths shorter than eight times the minimum resolvable 
wave and slightly degraded fit to the longer wavelengths. 
Therefore, in the application of the Barnes method to irregularly 
spaced data, successively decreasing the influence of the more 
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distant observations is still advisable if longer wavelengths are 
present in the field of data. 

It also was found that no single selection of filter 
parameters for the two-pass Barnes method gives the best analysis 
for all wavelengths. A three-pass hybrid method is shown to reduce 
this problem. 
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Abstract 

The use of objectively analyzed fields of meteorological data 
for complex diagnostic studies and for the initialization of 
numerical prediction models places the requirements upon the 
objective method that derivatives of the gridded fields be accurate 
and free from interpolation error. A modification of an objective 
analysis developed by Barnes provides improvements in analyses of 
both the field and its derivatives. Theoretical comparisons, 
comparisons between analyses of analytical monochromatic waves, and 
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comparisons between analyses of actual weather data are used to 
show the potential of the new method. The new method restores more 
of the amplitudes of desired wavelengths while simultaneously 
filtering more of the amplitudes of undesired wavelengths. These 
results also hold for the first and second derivatives calculated 
from the gridded fields. Greatest improvements were for the 
Laplacians of the height field; the new method reduced the variance 
of undesirable very short wavelengths by 72 percent. Other 
improvements were found in the divergence of the gridded wind field 
and near the boundaries of the field of data. 
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